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ABSTRACT 

Existing  experimental  results  of  hypervelocity  impact  tests  have  been 
gathered  from  various  sources  and  the  composite  data  are  presented  and 
discussed.  The  extrapolated  results  to  higher  velocities  are  seen  to  differ 
with  the  theoretical  prediction  from  the  perfect  fluid  model  proposed  by 
R,  L,  Bjork,  Re-examination  of  his  assumptions  and  the  experimental  results 
indicate  the  desirability  of  including  the  effects  of  the  viscoS’ty  and  the  yield 
strength  of  the  materials  into  the  mathematical  model.  A  visco-plastic  model 
for  hypervelocity  impact  is  then  formulated  to  meet  these  requirements.  This 
is  accomplished  by  introducing  a  viscosity  factor  (iQ  and  a  yield  stress  TQ 
into  the  perfect  fluid  equations. 

The  equations  governing  the  visco-plastic  model  are  then  studied  and 
the  characteristic  features  of  the  theory  are  deduced.  Certain  dimension¬ 
less  parameters  are  found  which  determine  the  relative  importance  of  the 
inertial,  viscous  and  plastic  effects  during  the  various  stages  of  the  hyper¬ 
velocity  cratering  process. 

To  exhibit  quantitatively  the  importance  of  including  the  viscous  and 
plastic  effects,  a  one-dimensional  impact  model  was  studied.  In  this  study 
the  values  of  and  Tq  are  varied  since  definitive  data  are  available  for 
neither  parameter  in  the  hypervelocity  range.  Two  distinct  finite  difference 
schemes  have  been  developed  for  performing  the  required  calculations  on  an 
IBM  7090  computer.  These  are  described  in  detail. 

The  results  of  the  calculations  are  related  to  the  qualitative  model  of 
crater  formations  that  has  evolved  from  experimental  studies  in  which  the 
actual  cratering  process  has  been  monitored.  It  is  concluded  that  the  viscous 
and  strength  effects  strongly  affect  the  cavitation  process  which  is  the 
essential  mechanism  of  crater  formation.  Finally,  experiments  are  suggested 
which  would  provide  the  necessary  data  to  verify  and  extend  the  theory. 
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INTRODUCTION 


For  satellites  and  for  manned  space  flight,  the  presence  of  meteoritic 
material  in  space  accounts  for  one  of  the  important  environmental  hazards. 

The  velocity  of  the  meteoroids  is  estimated  to  range  from  0.  5  to  7.  2  cm/ 
microsecond.  *  A  similar  problem,  potentially,  is  the  damage  that  may  be 
inflicted  on  a  ballistic  missile  by  high  speed  particles  produced  and  directed 
by  artificial  means.  For  these  reasons,  increasing  attention  has  been  re¬ 
cently  focused  on  the  problems  attendant  to  the  collision  of  a  projectile  and  a 
target  in  the  hypervelocity  regime.  The  ultimate  objective  of  these  studies 
is  to  determine  the  minimum  hull  thickness  required  to  ensure  that  the  space 
vehicle  is  not  pierced.  More  explicitly,  the  dependence  of  the  minimum  thick¬ 
ness  on  the  various  parameters  is  the  information  sought. 

Until  quite  recently,  laboratory  techniques  were  only  capable  of  pro¬ 
ducing  velocities  of  less  than  about  0.  6  cm /microsecond.  A  great  deal  of  data 
are  available  describing  the  craters  formed  at  these  lower  speeds.  In  the 
past  year  a  method  has  been  developed  for  projecting  hypervelocity  pellets  up 
to  2  cm /microsecond.  At  present,  however,  little  data  are  available  in  this 
range.  As  an  adjunct  to  such  experiments,  empirical  formulas  and  cratering 
theories  have  been  proposed  to  extrapolate  the  experimental  results  to  cover 
the  velocity  range  of  interest.  In  the  following  both  the  data  and  the  theories 
are  reviewed  and  the  conclusions  used  as  a  basis  for  the  formulation  of  a 
visco-plastic  model  for  hypervelocity  impact  which  takes  into  account  the 
viscosity  and  strength  of  the  projectile  and  target  materials  as  well  as  their 
com  press  ibilitie  s . 

The  visco-plastic  model  is  proposed  only  for  materials  that  behave  in  a 
ductile  manner  when  impacted.  The  phenomenon  is  quite  different  for  brittle 
materials  such  as  rock.  The  survey  of  the  experimental  data  is  restricted  to 
thick  metal  targets  impacted  at  normal  incidence.  The  use  of  a  target  which 
is  essentially  a  semi-infinite  body  rules  out  geometrical  complications  such 
as  reflected  shockwaves  at  free  surfaces,  bending,  etc.  The  large  majority 
of  past  investigations  have  been  similarly  restricted. 


The  gram -centimeter-microsecond  system  of  units  is  used  throughout  this 
report. 
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REVIEW  OF  EXISTING  DATA  AND  THEORIES 


When  velocity  limits  of  power  propellant  guns  were  reached 
experimenters  provided  the  first  hypervelocity  impact  data  through 
the  utilization  of  light  gas  (hydrogen,  helium)  projectors.  The  first 
light  gas  projector  was  developed  and  constructed  at  the  New  Mexico 
School  of  Mines  by  Dr.  W.  D.  Crozier  and  Dr.  William  Hume.  Most 
of  the  hypervelocity  data  presently  available  have  been  obtained  by 
use  of  light  gas  guns.  The  upper  limit  reported  as  attained  by  them  is 
about  0.  6  cm/microseconds. 

(a)  Survey  of  Data 

The  data  which  will  be  presented  have  been  abstracted  from  the 
reports  of  tests  using  both  these  high  velocity,  single  particle  impact 
techniques,  References  1  through  14.  The  particle  properties  are 
therefore  well  defined.  Our  object  is  to  analyze  the  data  for  charac¬ 
teristic  features  of  the  impact  phenomenon;  results  are  not  to  be 
included  if  the  particular  projectile  and  target  material  combination 
has  been  studied  in  a  single  series  of  tests. 

Changes  in  the  crater  profile  as  the  velocity  is  increased  over 
the  known  range  are  illustrated  by  Figures  1  through  4.  The  ratio  of 
the  penetration  to  the  diameter  of  the  crater,  Pc/Dc,  is  plotted  for 
various  projectile  materials  impacting  massive  targets  of  lead,  copper, 
steel,  and  aluminum  respectively.  The  low  velocity  scatter  is  assoc¬ 
iated  with  undeformed  projectile  penetration.  At  higher  velocities, 
depending  on  the  strength  of  the  projectile  and  target,  the  projectile 
deforms  plastically.  As  the  impact  velocity  increases  still  further, 
these  data  show  that  for  most  projectile  materials  the  crater  profile 
parameter,  Pc/Dc,  approaches  0.5,  the  value  corresponding  to  hemi¬ 
spherical  craters. 

For  some  material  combinations  the  velocity  at  which  the  0?  5  level 
is  attained,  if  attained  at  all,  is  seen  to  be  quite  high,  especially  for 
cases  in  which  the  yield  strengths  of  the  projectile  and  target  are  high. 
This  is  demonstrated  by  Figure  3,  which  shows  that  tungsten  carbide 
and  aluminum  alloy  projectiles  impacting  steel  have  not  Attained  Pc/Dc 
0.  5  for  velocities  up  to  0.  5  cm/ microsecond.  In  the  former  the  ratio 
is  closer  to  about  0.  6,  while  it  is  closer  to  0.  4  in  the  latter. 


2 


In  Figures  5  through  8  the  penetration  parameter,  Pc/Dg,  has  been 

plotted  as  a  function  of  the  impact  velocity  for  results  where  projectile  and 
target  materials  were  the  same.  Here  Ds  denotes  the  diameter  of  a  sphere 
having  the  same  mass  as  the  actual  projectile.  Pc  again  denotes  the  depth 

of  a  crater.  The  impact  ve!  Deity  is  expressed  in  dimensionless  terms  through 
division  by  the  velocity  of  sound  in  the  undisturbed  target  material,  v0/cQ. 

When  both  projectile  and  target  are  of  low  strength  the  log-log  plot  of  the  data 
is  seen  to  be  well  fitted  by  a  straight  line  of  slope  2/3,  Figures  7  and  8.  For 
materials  of  greater  strength  the  slope  at  low  velocities  differs  but  appears 
to  tend  towards  2/3  at  the  higher  velocities,  Figures  5  and  6. 

In  Figures  9  through  12  the  crater  volume  parameter,  Vc/Vp,  i8  plotted 

against  the  dimensionless  velocity  for  the  projectile  and  target  materials 
corresponding  to  Figures  5  through  8  respectively.  Here  Vc  and  Vp  denote 
the  crater  and  projectile  volumes  respectively.  A  straight  line  of  slop  2  is 
seen  to  be  in  excellent  agreement  with  the  plot  on  log-log  paper  in  all  cases, 
even  at  the  lower  velocities.  This  result  is  compatible  with  the  variation  of 
Pc/Ds  with  (v0/c0 )2/3  and  with  the  tendency  of  the  craters  to  approach  a 
hemispherical  shape  as  described  above. 

(b)  Empirical  Formulas 


Nearly  everyone  who  has  obtained  experimental  results  has  developed 
empirical  formulas  that  represent  his  observations  to  a  fair  degree  of 
accuracy.  The  formulas  are  not  based  on  rational  theoretical  grounds  and 
generally  are  not  valid  for  velocities  and  projectile-target  material  com¬ 
binations  outside  the  restricted  regions  covered  by  the  tests.  Recently, 
Bruce  (Ref.  15)  has  fitted  formulas  to  the  composite  firing  results  which 
are  consistent  with  the -assumptions  that  in  the  hypervelocity  range  Pc/Dc= 

0.5,  Pc/Dg  ~(v0/c0)2^3  and  Vc /Vp  (v0/c0)2.  His  are  probably  one  of  the 
two  moat  reliable  sets  of  engineering  formulas  presently  available,  but  even 
then  extrapolation  to  other  materials  or  to  velocities  greater  than  about  1. 0 
cm /microsecond  is  hazardous. 

More  recently,  Herrmann  and  Jones  (Ref.  16)  have  also  gathered  and 
analyzed  the  available  data  on  cratering  by  high  speed  impact  in  semi¬ 
infinite  targets.  They  found  that  the  penetration  in  the  high  velocity  region 
is  best  fitted  by  the  dimensionless  equation 

(1)  Pc/Ds=  k  (Pp/Pt)2/3  H  where  H-  Pt  vQ2/ht. 

Here  k  is  a  constant  near  0.  36  for  most  materials;  Pp  and  P^are  the  densities 
of  the  projectile  and  target  materials;  1^  is  the  Brinnell  Hardness  of  the 
target.  The  Brinnell  Hardness  is  defined  as  the  load  applied  to  a  spherical 
stylus  divided  by  the  area  of  the  resultant  depression,  and  thus  has  the 
dimensions  of  a  stress. 
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There  is  an  approximate  dependence  of  indentation  hardness  on  the 
shear  yield  stress,  rQ.  The  relation  is  approximately 

(2)  h»7.2T0. 

It  would  therefore  appear  that  the  strength  of  the  target  material  is  an  im¬ 
portant  factor  in  the  cratering  mechanism.  In  fact  Herrmann  and  Jones  tried 
various  parameter  combinations  in  which  the  material  strength  was 
neglected,  and  found  that  none  could  fit  the  composite  experimental  data. 

(c)  Previous  Cratering  Theories 

A  number  of  simple  theories  of  hypervelocity  impact  have  been  proposed 
for  the  prediction  of  penetration  (Refs.  17  through  23).  Most  of  the  theories 
show  little  agreement  with  the  experimental  data;  this  is  not  surprising  since 
neither  realistic  compressibilities  nor  flow  geometries  are  introduced  into 
any  of  the  derivations. 

The  only  serious  attempt  to  calculate  the  phenomenology  of  hypervelocity 
impact  from  basic  physical  equations  has  been  made  by  Bjork  (Ref.  24).  He 
uses  a  hydrodynamical  model  and  treats  the  rotationally  symmetric  case  of 
two  dimensional,  unsteady,  compressible  flow  in  a  semi-infinite  target 
under  normal  impact  by  a  cylindrical  projectile  of  the  same  material.  In 
setting  up  the  mathematical  model  Bjork  assumes  that  a)  the  elastic  waves 
can  be  neglected  since  the  stresses  and  particle  velocities  carried  by  these 
waves  are  much  less  than  those  caused  by  shock  waves,  b)  the  flow  is  strictly 
adiabatic,  c)  the  strength  of  the  target  and  projectile  material  is  negligible, 
and  d)  the  flow  is  inviscid. 

The  results  of  Bjork's  calculations  are  summarized  by  the  equations 

AlonAl:  Pc/D,  *  2.  09  (v0/c0)1/3 

(3)  ,,, 

FeonFe:  Pc/Ds=  1  69  (vQ/c o)1'  • 

Within  the  numerical  error  of  the  computational  procedure  he  found  the 
craters  to  be  hemispherical  and  so,  from  geometry  alone, 

(4)  Vc/Vp=  4  (Pc/D8)3. 

Substitution  of  equation  (3)  into  (4)  yields 

(5)  AlonAl  :  Vc/Vp  =  36.  5  v0/c0 

FeonFe  :  Vc/Vp  =  19.  3  v0/cQ. 

According  to  this  model,  therefore,  the  penetration  at  hypervelocities 
actually  increases  with  about  the  one-third  power  of  velocity  rather  than  with 
the  two-thirds  power  as  extrapolation  from  the  lower  velocity  tests  would 
indicate. 
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Since  densities  and  velocities  are  the  primary  variables  in  the  perfect 
fluid  model,  equations  (3)  and  (5)  should  apply  equally  well  to  other  materials 
of  equal  density.  The  predicted  curves  for  the  penetration  parameter  and 
crater  volume  parameter  are  superimposed  on  the  experimental  plots  of 
Figures  5,  6,  and  Figures  9,  10  respectively.  The  mismatch  even  at  the 
higher  velocities  in  the  experimental  range  causes  one  to  have  reservations 
in  using  the  model  to  predict  results  in  the  upper  reaches  of  the  0.  5  to  7.2 
cm /microsecond  range.  Therefore,  this  leads  to  examination  of  his  four 
basic  assumptions.  Since  a)  and  b)  are  almost  certainly  justified  attention  is 
focused  On  c)  and  d). 

(d)  Viscous  and  Strength  Effects 

As  to  neglecting  the  strength  of  the  material,  for  dynamic  conditions 
where  the  duration  of  loading  is  very  small  (say  in  the  order  of  microseconds) 
the  dynamic  yield  strength  may  rise  by  several  orders  of  magnitude.  Al¬ 
though  suggested  by  earlier  workers  (Refs.  25  and  26)  and  emphasized  above, 
this  dependence  of  the  cratering  process  on  the  material  strength  has  recently 
been  directly  illustrated  by  an  experiment  carried  out  at  the  Carnegie  Institute 
of  Technology  (Ref.  13).  The  experiment  consisted  simply  of  firing  steel 
pellets  (0.  18  gram,  0.  5  cm/microsecond)  into  targets  of  lead,  cadmium,  zinc, 
and  copper  and  observing  the  craters  produced  by  the  impact  as  the  target 
temperature  was  varied  over  a  wide  range.  The  crater  volume  plotted 
against  temperature  showed  abrupt  changes  at  certain  temperatures  that  are 
identified  with  similar  discontinuities  in  tensile  tests  made  on  the  metals  as 
the  temperature  was  varied.  The  critical  temperatures  are  the  points  where 
certain  metallurgical  changes  such  as  transition  from  brittle  to  ductile  be¬ 
havior  or  stress  anneal  occur. 

That  the  viscosity  of  the  target  plays  an  important  role  in  impact  phenom¬ 
enon  may  be  seen  from  post  mortem  metallurgical  examination  of  the  micro¬ 
structure  surrounding  a  high  velocity  crater  (Ref.  2).  Practically  no  change 
in  the  shape  of  the  grains  in  the  material  occurs  under  the  crater  where  the 
strain  rate  is  small.  On  the  other  hand,  the  grains  are  found  to  be  elongated 
considerably  along  the  sides  of  the  crater  where  the  strain  rate  is  greatest. 

The  distortion  is  caused  by  the  shear  stresses  which  result  from  the  high 
strain  rate.  Bjork's  inviscid  model  allows  only  for  hydrostatic  pressure,  and 
thus  does  not  take  this  effect  into  account.  The  viscosity  is  also  important 
from  the  point  of  view  that  without  considering  it  there  is  apparently  no  way 
of  introducing  anisotropic  stresses  into  the  flow.  This  must  be  done  if  a 
strength  effect,  which  is  clearly  essential,  is  to  be  introduced. 
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VISCO- PLASTIC  MODEL 


In  this  section  the  problem  will  be  considered  anew  in  order  to  formulate 
a  mathematical  model  which  takes  into  account  both  the  strength  and  viscosity 
of  the  materials  involved. 


(а)  Bingham  Model 

When  an  ultra -high- speed  projectile  strikes  a  target  a  strain-rate  which 
depends  on  the  impact  velocity  is  imposed  on  the  projectile  and  target 
materials.  The  plastic  deformation  will  be  resisted  not  only  by  the  static 
yield  stress,  but  also  by  viscosity  stresses  with  magnitudes  dependent  on  the 
strain  rate.  Both  these  effects  seem  to  be  important.  Such  a  plastic  solid 
(exhibiting  visco-plastic  flow)  is  most  simply  represented  by  a  Bingham 
model.  The  material  is  considered  rigid  if  stressed  below  its  yield  strength, 
whereas  above  this  value  the  material  acts  like  a  Newtonian  viscous  liquid; 
a  schematic  representation  of  such  a  material  is  given  in  Figure  13(a).  In 
simple  shearing  flow,  in  which  the  velocity  is  q  (y),  and  dq/dy  is  a  constant 
D,  this  means  that 

t-T0=Mo  D  (T£T0) 

(б)  t+t0=m0d  r*-r0) 

D=0  (  |T|<T0>, 

where  is  the  yield  value  of  the  shearing  stress  and  MQ  is  a  constant  (if 

temperature  dependence  and  pressure  gradient  effects  are  neglected)  deter¬ 
mining  the  magnitude  of  the  strain-rate  effect. 

Now  the  usual  definition  of  the  viscosity  of  a  liquid  is  based  on  Newton's 
assumption  that  the  shear  stress  and  the  strain- rate  are  related  according  to 

T=  M  D, 

where  the  viscosity  coefficient  M  is  a  constant  at  constant  temperature.  When 
equations  (6)  are  written  in  this  form,  the  result  is 


(7)  T=  M( D)  D, 


where  the  strain-rate  dependent  viscosity  coefficient  M  =  M(D)  is  given  by 


(8)  M(D)  =  M0 

=© 


(W>  Tq) 
(Kl<  T0). 


This  dependence  of  M  on  D  is  illustrated  in  Figure  13(b).  The  viscosity 
coefficient  decreases  with  increased  strain-rate. 
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The  use  of  the  Bingham  model  is  proposed  to  bridge  the  transition  from 
the  plastic  to  the  hydrodynamic  regimes  (Ref.  27).  It  is  analogous  to  Malvern's 
model  which  seems  to  have  successfully  bridged  the  transition  from  the  elastic 
to  the  plastic  regime  (Ref.  28).  Malvern  assumed  that  the  strain-rate  was 
directly  proportional  to  the  difference  between  the  instantaneous  stress  and 
the  static  stress  corresponding  to  the  strain,  so  that  for  longitudinal  stress 
in  a  slender  rod  the  relation  takes  the  form 

(9)  -^-[EC-  cr]=£[a  .  F  (e)] 

where  E  is  the  elastic  modulus;  O  and  f  are  the  instantaneous  values  of 
longitudinal  stress  and  strain;  and  F  (C)  defines  the  static  stress -strain  curve. 
These  quantities  are  depicted  in  Figure  14.  As  €  increases  ^becomes  negligible 
compared  to  E€,  and  F  (f)  levels  out  and  may  be  approximated  by  a  constant 
Cq.  When  very  large  strains  are  involved  eauation  (9)  may  therefore  be 
approximated  by  the  relation 

(10)  a  -  CTq  ~v(  (<7 £<7o) 

The  similarity  of  equations  (6)  and  (10)  is  obvious. 

(b)  Axisymmetric  Stress  to  Strain-Rate  Relation 


In  order  to  extend  the  idealized  visco-plastic  model  to  a  material  sub¬ 
jected  to  hypervelocity  impact  the  stress  vs.  strain- rate  relations  (7)  and  (8) 
must  be  generalized.  As  formulated  there,  only  the  properties  of  the  material 
in  simple  shear  flow  are  defined;  corresponding  relations  for  three  dimen¬ 
sional  flow  are  required.  The  expressions  will  be  somewhat  simplified  by 
assuming  that  the  projectile  is  axisymmetric  in  geometry,  and  that  it  strikes 
a  semi-infinite  target  normal  to  the  axis  of  symmetry.  Furthermore,  the 
angular  momentum  of  the  projectile  is  assumed  to  be  zero  at  impact. 

Under  these  assumptions  the  only  non-vanishing  components  of  the  stress 
tensor  are  the  normal  stresses  rTI,  TS6  •  Tzz  anc*  the  shearing  stress  Trz  = 

TZr.  The  tensor  may  be  decomposed  into  two  parts,  a  component  pro¬ 
ducing  distortion  and  a  component  producing  only  a  volume  change: 

<">Tik“  <Tik *  eaik>  -  f>4ik 
•  Tik  -  P8ik 

where  6^  is  the  Kronecker  delta  and 

(12)  ’Tk*  Tik  +  P5ik;  P=-jTii- 

* 

rik  is  associated  with  the  internal  friction  of  the  medium  and  p  is  assumed 
equal  to  the  hydrodynamic  pressure.* 


*This  assumption  is  equivalent  to  setting  the  "second  coefficient  of  viscosity" 
M'  *  -2M  /3 .  See  Ref.  29. 
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Flow  is  assumed  to  occur  whenever  the  von  Misea  criterion  is  satisfied, 
i.e.  whenever 


(13)  j  2  ^  t2 
o 


where  T  = 


r*2  +  -L  (t*2  +  t*2  +  r *2)  . 

rz  2  rr  gg  zz' 


T  is  a  measure  of  the  magnitude  of  the  components 

* 


tensor,  T  . 

lk 


of  the  distortion  stress 


The  non-vanishing  components  of  the  strain-rate  tensor  are  defined  by 

(14)  oq  dq 

D  =  2  -nrX  D  =  —  q  D  =2  -?-*■ 
rr  °r  eg  r  ^r  zz  oz 

dqrx  b% 

D  =  D  =  r— r  +  r-2- 
rz  zr  d  z  dr 


This  tensor  may  be  decomposed  into  a  distortional  strain-rate  and  a  dilational 
strain-rate: 


Dik  = 


D. 


ik 


2/1  dp\ 

3 \P  dt  J 


ik 


Here-(1/P)  dP/dt  is  the  volume  strain  rate,  the  factor  1/3  is  required  to  obtain 
the  linear  strain  rate,  and  the  2  enters  because  of  our  definition  of  D..  The 
continuity  equation  for  compressible  flow,  11 


(15)^+  P  div  q  = 


0. 


may  be  used  to  write  the  above  equation  in  the  form 

<161  Dik-  D‘k  +  ldi^  6ik- 

To  generalize  (7),  the  distortional  stress  and  strain-rate  components  are 
related  by 


<I7>rU-',DU 

(tZ^tZ) 

O 

$ 

2  2 

D.,  =  0 
ik 

(T  <  V 

where  the  dependence  of  P  on  the  distortional  rate  of  strain  must  generalize 
(8)  in  a  natural  manner.  Substituting  (12)  and  (16)  into  (17)  gives 
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(18) 


Tik  "  "p5ik  +  M  °ik  '  T  4  V 

If  Mwere  a  constant,  equation  (18)  would  reduce  to  the  classical  stress  to 
strain-rate  relationship  of  the  Navier-Stokes  theory  of  hydrodynamics.  On 
the  other  hand,  if  one  does  not  assume  that  the  mean  normal  stress  is  the 
thermodynamic  pressure,  (12),  the  relationship  is  complicated  by  the  factor 
-2P/3  in  the  last  term  being  replaced  by  the  second  coefficient  of  viscosity Ji1  • 


Since  the  viscosity  coefficient  p  depends  on  the  distortional  strain-rate  it 

* 

must  be  a  scalar  function  of  the  three  invariants  of  the  tensor  D  i.e.  M 

must  be  independent  of  the  particular  frame  of  reference.  Now  the  first 
invariant  is  zero  and  the  dependence  on  the  third  invariant  is  small.  There¬ 
fore,  the  viscosity  may  be  assumed  to  depend  only  on  the  second  invariant: 


(19) 


d2  = 


>*2 

rz 


+  1 


<D*2  + 
rr 


*2 

’ee 


+  D*2) 
zz' 


=  D2  +1 
rz  2 


(D2  + 
r  r 


Dee+  dmi  -  \  |div  4*2- 


Thus,  the  natural  generalization  of  the  Bingham  model  to  the  case  of 
axisymmetric  flow  is  for  the  viscosity  M  to  depend  on  the  rate  of  strain 
according  to  functional  relation  (8)  with  D  given  by  (19): 


(20) 


M  =  + 


M=  CO 


+  ?  (°2  +  °2 


rr 


98  +  DL>  *  7  <div 


iv  q)2J 


172 


r^r2) 


(r2<  r2). 


The  generalization  of  the  Bingham  model  to  3  dimensions  has  been  formulated 
by  Oldroyd  (Ref.  30)  for  the  case  of  rectangular  coordinates. 

(c)  Formulation  of  Governing  Equations 

Now  that  we  have  the  expression  for  the  dependence  of  the  viscosity  on 
the  strain- rate  the  principle  of  the  conservation  of  momentum  may  be  applied 
to  write  the  two  equations  of  motion  corresponding  to  the  Navier-  Stokes 
equations  for  axially  symmetric  flow.  The  result  is  the  same  except  for  the 
non-constancy  of  M.  These  two  equations  together  with  the  continuity  equation 
(15),  which  results  from  the  conservation  of  matter,  give  three  equations  in 
the  four  dependent  variables  q  ,  q  ,  p,  P.  To  write  another  equation  some 

assumption  that  specifies  the  particular  type  of  flow  must  be  introduced.  An 
assumption  is  made  identical  to  that  of  Bjork,  that  no  heat  is  transferred 
between  neighboring  particles  of  the  material.  This  leads  to  the  well  known 
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energy  equation  of  hydrodynamic®  (Pef.  31).  This  equation,  however,  in¬ 
troduces  a  fifth  dependent  variable,  the  internal  energy  per  unit  mass  U. 

This  difficulty  is  circumvented  by  assuming  that  an  experimentally  determined 
equation  of  state  relates  this  thermodynamic  parameter  to  the  other  state 
parameters  p,  P.  A  discussion  of  our  choice  for  the  equation  of  state  is 
given  later  in  the  report. 

The  final  five  equations  are  as  follows: 


(21) 


(22) 


(23) 


(24) 

(25) 


(Mass) 


(Radial 

Momentum) 


(Axial 

Momentum) 


P,  UP.  u) 


P 


dU+  d  (1/P) 
dt  P  dt 


P  (D)  D2. 


(State ) 
(Energy) 


The  material  time  derivatives  may  be  expanded  according  to  the  relation 


+  q  •  grad  ( 


). 


and  the  divergence  for  the  axially  symmetric  case  is 


1  d  d 

div  q  =  -  —  (rq  )+  — —  q  ■ 
^  r  Jr  '  3z  z 
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STUDY  OF  GOVERNING  EQUATIONS 


The  introduction  of  the  effects  of  viscosity  and  strength  into  the  perfect 
fluid  equations  has  necessarily  led  to  complication  of  the  partial  differential 
equations  governing  the  flow.  As  the  next  logical  step  these  equations  will 
be  studied  qualitatively  to  determine  the  dimeneionless  parameters  which 
control  the  relative  importance  of  the  inertial,  viscous  and  strength  effects 

(a)  Characteristic  Numbers 

Consider  two  distinct  projectile- tar  get  systems,  both  having  projectile 
and  target  geometrically  similar .  Choose  the  diameter,  D#,  of  the  sphere 

with  mass  equal  to  that  of  the  projectile  as  the  characteristic  length-  We  may 
set 

<26>  Ds2=KlDsr 

The  same  relation  holds  between  all  other  pairs  of  corresponding  points  if 
cratering  flows  are  dynamically  similar: 

r2  =  Ki  rl  zZ~  K1  V 

For  two  corresponding  times  t^  and  t^  set 

(27)  t2  =  K2tr 

Since  velocity  means  traveling  a  certain  length  in  a  certain  time,  the  charac¬ 
teristic  velocities  (impact  values)  in  the  two  systems  are  related  according  to 

(28)  Kx 

V2  =  *1  "1 


and  hence 


K 


qr2  ~  K,  qrl 


*1  o 

qz2  =  kT  q 


zl. 


Furthermore,  we  may  set 


(29) 


P2  =  k3  Pj 


*02=  K6*ol 


p2  "  K4  p  1 


To2=K7Tol- 


U2  =K5  U1 


We  now  use  relations  (26)  through  (29)  to  express  the  magnitudes  of  the 
variables  of  the  second  cratering  flow  in  terms  of  those  of  the  first.  This  is 
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done  in  each  of  the  equations  governing  the  visco-plastic  formulation.  The 
K  factors  for  the  two  terms  in  (21)  are  the  same  so  that  the  continuity  equations 
for  the  two  flows  are  identical.  The  viscosity  M  =  M(D)  which  appears  in  the 
momentum  and  energy  equations  is  homogeneous  in  K  only  if 

K2 

(30)  K6=K2K?  or  K6=K1K?|^. 

Then  M2  =  M  ^  and  all  terms  in  (22)  and  (23)  have  the  same  factor  provided 


(31) 


K1K3 


K- 


K6  , 

K1  *2  Kl‘ 


The  specific  energy  term  (consisting  of  kinetic  and  internal  energy  com¬ 
ponents)  in  (25)  is  homogeneous  in  K  only  if  Kg  =  K^/K^  or,  from  (31), 


(32) 


K1 

K  =  — K 


K„ 


Then  all  the  terms  in  the  energy  equation  have  the  same  factor  provided 
(33)  K3K5  K, 

K2  k: 


c2- 

c 


Upon  substitution  from  (26),  (27),  and  (29)  the  required  K  relations  may 
be  rewritten  in  dimensionless  forms.  Relation  (30)  is  equivalent  to 


Dsl  Tol  _  Ds2  To2 


Mol  V! 


Mo2  v2 


Relations  (31)  are  equivalent  to 
D 


si 


1 


Ds2  P2  v2 


"ol  ^02 

Relation  (32)  is  equivalent  to 

Ui_  Ug_ 

2  "  2 
V1  v2 

and  relation  (33)  gives  nothing  new. 


and 


P  v  ‘ 

*  1  1 


P2  w2 
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Necessary  conditions  for  dynamic  similarity  of  the  two  systems  are  that 
they  have  the  same  values  for  the  dimensionless  ratios,  Reference  32 


(34) 


B 


r 

o 

V 

o 


D  p 

s  o 


V 

o 


Here  the  subscript  o  denotes  that  the  quantity  is  evaluated  in  the  un¬ 
disturbed  state.  Bq  and  Rq  are  the  familiar  Bingham  and  Reynolds  numbers 

which  ordinarily  arise  separately  in  the  theory  of  slow  visco-plastic  flow  and 
the  theory  of  viscous  liquids  respectively.  Here  both  occur  as  the  model  in¬ 
cludes  the  two  types  of  flow.  M*  is  the  generalized  Mach  number  which  for 
an  ideal  gas  may  be  decomposed  into  two  factors,  the  ratio  of  specific  heats 
and  the  ordinary  Mach  number  M.  It  is  interesting  to  note  that  if  the  charac¬ 
teristic  pressure  pQ  is  taken  to  be  the  Brinnell  Hardness,  see  equation  (1), 

then 


m*=>/h7 


The  characteristic  specific  internal  energy,  U^,  may  be  taken  to  be  the  energy 
required  to  melt  or  vaporize  a  unit  mass  of  the  medium. 

Since  the  governing  equations  also  include  the  equation  of  state,  equality 
of  numbers  (34)  are  not  sufficient  to  ensure  that  two  geometrically  similar 
flows  are  dynamically  similar.  The  equation  of  state  does  not  lend  itself  to 
this  type  of  investigation  since  it  is  an  empirical  formula. 

(b)  Relation  to  Perfect  Fluid  Model 


The  equations  governing  the  flow  can  be  written  in  dimensionless  form  by 
setting 


(35) 


r  *  Lr ' 

z  =  Lz' 

p=  p  P' 
o 

q  -  vq  ' 
ir  ^r 

%  -vV 

U  =v2U' 

.  2  , 
p  =PqV  p' 

t  =Lv'*t' 

R  -LP  vM*1 
o  o 

M  -LP  vR'1 
o  o 

T  =p  v2BR'1 

O  O 

B  =  L  r  (M  V 

0  O  I 
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where  the  dimensionless  quantities  are  indicated  by  a  prime.  The  governing 
equations  then  become 


(36) 


dp' 

dt 


+  P'  div‘  q'  =  0 
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The  left  side  of  (40)  may  also  be  written  as 


upon  substitution  from  (36). 

The  quantities  in  the  brackets  in  each  of  equations  (37),  (38),  (40)  will  be 
of  order  unity  at  each  point  in  space  and  time  that  L,  v  are  truly  representative 
values.  If  so,  and  if  one  of  the  three  numbers  R,  1,  B  is  much  greater  than 
the  other  two,  a  single  set  of  terms  in  each  equation  predominates,  and  the 
motion  may  be  approximately  described  by  equating  the  coefficient  of  the  pre¬ 
dominant  number  to  zero  For  example,  if  R^l,  B,  the  equations  (37),  (38), 
(40)  reduce  to 


dq*r  dn. 

O' - -  +  — i— 

M  dt'  Sr' 


0 


dq'  a  . 

i - -  +  — tL.  s  o 

dt'  3z' 


+  P'  divq'«  0, 

respectively. 

Actually,  these  combined  with  (36)  and  (39)  are  the  perfect  fluid  equations, 
and  are  seen  to  be  valid  when  the  inequalities  RJM,  B  hold,  i.  e.  only  when 
the  inertial  effect  predominates.  This  will  be  the  case  during  those  parts  of 
the  cratering  process  when  the  characteristic  velocity  v  appearing  in  R  and 
B  is  high  enough,  i.  e.  ,  by  (35), 

(42)  v  » maximum  {y^QtP Q>  MQ/L 

As  the  rate  of  flow  decreases  so  does  the  characteristic  velocity  until 
eventually,  when  the  flow  ceases,  v  -  0.  Clearly,  in  the  late  stages  of  flow 
the  inequalities  required  for  the  perfect  fluid  approximation  to  be  valid  are 
reversed,  i.e.  ,  R«l,  B.  At  other  stages,  two  or  possibly  all  three  of  the 
groups  of  term 8  may  be  equally  important. 

It  is  seen  that  the  strength  terms  are  insignificant  during  the  early  stages 
of  flow  but  become  predominant  in  the  latter  stages.  The  question  arises  as 
to  the  importance  of  the  viscous  terms.  Clearly,  they  would  always  be 
negligible  only  if 

(43)  maximum  (R,  B)»l 

for  the  entire  range  of  v,  (0,  vq).  Since  R  and  B  are  respectively,  monotone 


15 


increasing  and  decreasing  functions  of  v  the  left  side  of  (43)  is  attained  when 
R=  B,  i.  e.  ,  when  the  flow  velocity  is  given  by 

v*»yto/?o  • 

The  required  inequality  for  the  viscous  terms  to  be  negligible  throughout  the 
cratering  process  therefore  becomes 

<44>  • 

At  the  instant  of  impact  an  exceedingly  large  strain-rate  gradient  is  im¬ 
posed  at  the  surface  of  contact.  The  use  for  example,  of  the  diameter  of  the 
projectile,  Dg,  for  the  representative  length  L  cannot  be  justified;  a  much 

smaller  value  is  actually  required.  Accordingly,  the  inequalities  (42)  and  (44) 
are  very  restrictive.  The  viscosity  is  probably  never  really  negligible  in  the 
time  interval  immediately  after  impact. 

(c)  Choice  of  Parameter  Values 

The  principal  difficulty  in  testing  these  conclusions,  and  a  real  difficulty 
in  the  application  of  the  governing  equations  themselves,  is  in  assuming  the 
viscosity  coefficients  for  structural  materials  such  as  steel,  copper,  or 
aluminum.  The  viscosity  involved  here  has  little  relation  to  ordinary  creep 
data.  The  choice  of  an  appropriate  value  for  the  dynamic  shear  yield  strength 
T  is  nearly  as  troublesome.  It  is  known  to  increase  with  increasing  strain- 

rate,  but  values  are  not  known  in  the  hypervelocity  impact  range.  Therefore, 
to  check  the  philosophy  used  in  the  construction  of  the  mathematical  model 
exploratory  calculations  had  to  be  performed  in  which  the  values  of  and  T ° 

were  varied  over  several  orders  of  magnitude  and  the  impact  velocity  vq  was 

allowed  to  assume  a  set  of  values  in  the  range  of  interest. 

As  a  first  step  in  choosing  trial  values  for  MQ  the  visco-plastic  formula¬ 
tion  will  be  applied  to  the  case  of  a  long  rod  subjected  to  uniaxial  loading 
(along  z  -  axis).  Then 

D  =  0  D  =  Do o  =  -  4*D  div  q  =  0 
rz  rr  oo  2  zz  ^ 


and  the  viscosity  coefficient,  (20),  reduces  to 

r 


M  =  M  +  - ~~ 

°  V3|aqs 

■35" 
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Hence,  by  (12),  the  uniaxial  stress  and  strain- rate  are  related  according  to 


r 

ZZ 


3  M 


o  dz 


This  relation  compares  with  (10)  upon  setting 


(45)  n  /3. 

r 

Now,  Malvern  (Ref.  28)  has  applied  his  model,  (9),  t.o  experimental 
data  obtained  upon  subjecting  a  long  rod  to  a  plastic  wave.  In  this  way  values 
for  t|  in  the  elastic-plastic  transition  region  were  determined.  His  values 
may  be  substituted  into  (45)  to  obtain  approximate  values  of  for  the  plastic- 
hydrodynamic  transition  region: 


Steel:  M  =08  gm  cm"*  microsecond  ’ * 
o  ° 

Copper:  /iQ  =  0.4  gm  cm"'  microsecond'*. 


Others  have  also  obtained  approximately  the  same  values,  Reference  33.  Our 
calculations  have  all  been  for  iron  and  M  is  assumed  to  be  within  a  factor  of 
ten  of  0.8. 

For  mild  steel,  Reiner  (Ref.  34)  quotes  the  static  yield  stress  to  be 
approximately 


r 

o 


gm  cm 


-1 


microsecond’  , 


i.  e.  ten  kilobars.  For  our  calculations  the  dynamic  yield  stress  is  pertinent. 
The  above  value  is  assumed  to  be  the  lowest  value  likely,  and  is  varied  up 
to  one  hundred  times  as  large,  i.  e.  one  megabar. 

The  various  combinations  of  assumptions  for  Hq,  and  impact  velocity, 
v  ,  are  displayed  in  Table  I.  The  choices  of  vQ  (=0.  5,  4,  7.  5  cm /microsecond) 

represent  the  extremes  and  mean  of  the  meteoroid  velocity  range.  It  will  be 

desirable  to  have  further  results  for  v  =1.0  cm /microsecond.  The  values  of 

o 

the  dimensionless  parameters  Bq  and  Rq  are  also  listed  for  each  combination 

of  n  T  ,  v  .  They  are  seen  to  vary  widely  with  the  choice  of  combinations, 
o  o  o 
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TABLE  I 
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ONE  -  DIMENSIONAL  PROBLEM 


For  the  exploratory  calculations  we  have  considered  the  one -dimensional 
impact  of  two  semi-infinite  bodies.  Figure  15,  and  compared  the  results 
obtained  with  those  given  by  the  perfect  fluid  equations.  This  is  a  problem  of 
"plane  strain"  and  the  problem  of  a  long  rod  considered  above  is  a  problem  of 
"plane  stress". 

(a)  The  Eulerian  Formulation 


The  visco-plastic  equations  for  one -dimensional  flow  can  be  deduced  from 
the  axisymmetric  equations  by  equating  to  zero  qr  and  all  derivatives  with 

respect  to  r.  Then,  setting  q  as  q,  the  viscosity  coefficient  reduces  to 

z 


(46) 


4  =  4  + 
o 


n  =  oo 


) 


(T2  <  T 


2 

o 


). 


(47) 


The  components  of  the  distortional  stress  tensor 

*  -  *  2  ?q  _  *  4  ..6q 

i  =  -  -j  4  ;rr  T =T^dr 


Trr  T00 


bz 


zz 


T^az 


rz 


and  the  von  Mises  flow  statistic  therefore  reduces  to 


are,  by  (17), 

* 

-  0, 


(48) 


T 

o 


The  system  of  equations  governing  the  flow  reduces  to 


(51)  p  =  f  (P.  U) 


In  writing  (50)  we  have  used  the  relationship 
(53) 


To  show  that  this  is  the  case  is  equivalent  to  showing  that  the  velocity  gradient, 

dq/dz,  never  changes  from  just  above  zero  to  just  below  zero.  From  (46), 

(4?)  we  see  that  this  situation  cannot  occur  since,  for  example,  T  *  would 

_  zz 

have  to  change  from  just  above  -s/4/3  T  to  just  below  -  >/4/3  T  ,  involving  an 
impossible  discontinuity  of  stress.  °  ° 

(b)  The  Lagrangian  Formulation 

The  one  dimeneional  formulation  above  is  in  the  Eulerian  form.  The 
Lagrangian  form,  however,  gives  more  information  (it  keeps  track  of  each 
material  particle)  and  has  the  virtue  that  conservation  of  mass  is  automatic 
and  exact,  even  in  the  finite  difference  approximations.  Because  of  this 
inherent  greater  accuracy  the  Lagrangian  form  is  often  preferred  for  problems 
in  one  space  variable. 

In  Lagrangian  coordinates  the  motion  is  considered  in  terms  of  the  in¬ 
stantaneous  position  X  of  a  section  which  is  a  function  of  time  and  a  space 
coordinate  x  which  identifies  the  section.  Here  the  Lagrangian  coordinate  x 
will  be  taken  to  be  the  initial  position  of  the  section  at  time  zero.  The  velocity 
is  then  given  by 


The  transformation  from  the  Eulerian  variables  (z,t)  to  the  Lagrangian 
variables  (x,  t)  is  accomplished  by  the  relations 


(55) 

dg_/< 

ag  _p  /  9 

g\ 

>t/x 

az  po\S 

V/ 

where  g  is  any  dependent  variable.  In  writing  the  second  of  these  we  have 
used  the  fact  that  the  conservation  of  mass  implies 

SX_^o 

dx  p  • 
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In  Lagrangian  form  the  governing  equations,  corresponding  to  (49) 
through  (52)  respectively,  are  as  follows: 


(57)  > 


(58)  p  =  f  (P,U) 


(59) 


iu  1  m  _  4  P  /aq\*  /S'  To 

*  p*  W-J'V*  (-jy  +n  r 


o  ■  dx 


,6°)  g2 


The  viscosity  coefficient  becomes 

T 

q 

dx 

p  -  CD 

The  stress  components  become 


(T2>ro2) 


(t2<To2). 


(61) 


and 


(62) 


t  r  ^  ..  P  dq  T o  /dq\ 

Trr  =  T0e  =’P- P  a  "7F  Slgn  fa  ) 

o  dx  y/3  \dx/ 


rzz  =  -P  +T  ^  ^  +  »ign 


3  ®o  dx^3  °  *18”  (a*) 


T2  = 


=  1  u2  £_  /is\ 

J  P  £  \dx/ 


T  =0. 
rz 


o 


r 

o 


_p 

p 

o 
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(c)  Necessity  of  Numerical  Solution 


Analytical  methods  for  solving  the  governing  equations  are  not  available, 
even  for  these  one  dimensional  formulations.  It  was  therefore  necessary  to 
resort  to  finite  difference  techniques  for  the  computations.  To  cover  the 
desired  range  of  parameter  combinations,  listed  in  Table  I,  required  the 
development  of  two  separate  computational  schemes.  The  first  is  an  explicit 
difference  scheme  based  on  the  Lagrangian  formulation.  It  allows  the  two 
impacting  bodies  to  be  of  different  materials,  and  places  no  restriction  on 
For  large  values  of  and  v  ,  however,  it  is  advantageous  to  use  an 

alternate  implicit  scheme  based  on  the  Eulerian  formulation.  The  latter 
scheme  then  requires  less  machine  time,  but  it  is  only  valid  for  impact  be¬ 
tween  bodies  of  identical  material  with  M  >0.  The  two  methods  are  there¬ 
fore  mutually  complementary.  ° 

The  next  two  sections  describe  these  schemes  in  detail.  They 
may  be  omitted  by  the  more  casual  reader. 
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EXPLICIT  DIFFERENCE  SCHEME 


The  explicit  difference  scheme  represents  an  extension  of  a  method 
commonly  employed  for  perfect  fluids  to  allow  for  the  additional  terms  which 
occur  in  the  visco -plastic  model. 

(a)  Artificial  Viscosity 

In  a  perfect  fluid  shock  waves  occur  as  moving  surfaces  across  which  p, 

U,  p,  q  are  all  discontinuous.  The  differential  equations  (which  makes  no 
sense  on  such  surfaces)  must  then  be  augmented  by  jump  conditions  which 
serve  as  internal  moving  boundary  conditions  and  their  occurrence  vastly 
complicates  the  solution.  To  avoid  these  difficulties  von  Neumann  and 
Richtmeyer  (Reference  35)  introduced  a  purely  "artificial  viscosity"  which 
has  a  smoothing  effect  so  that  the  discontinuity  surface  is  replaced  by  a  thin 
transition  layer  in  which  the  dependent  variables  change  rapidly,  but  not  dis- 
continuously.  The  purely  artificial  dissipative  mechanism  was  chosen  of  such 
form  and  strength  that  the  required  smoothness  is  achieved  without  affecting 
the  flow  pattern. 

In  the  equations  which  are  being  investigated  here  a  true  viscosity  may 
or  may  not  be  present.  Such  a  smoothing  mechanism  must  be  included  if  all 
the  parameter  combinations  are  to  be  considered.  In  fact,  it  was  found  neces¬ 
sary  to  include  such  a  device  for  certain  cases  where  is  positive  but  small. 
The  real  viscosity  terms  are  proportional  to  the  strain-rate  whereas  the  arti¬ 
ficial  viscosity  (as  will  be  shown  in  (64)  )  is  quadratic  in  the  strain-rate.  The 
latter  is  therefore  more  effective  in  smoothing  extremely  strong  shocks.  On 
the  other  hand,  the  artificial  viscosity,  being  quadratic,  is  extremely  small 
in  the  smooth  part  of  the  flow  between  shocks,  where  one  wishes  the  true 
viscosity  coefficient  to  predominate. 


When  the  Lagrangian  equations  are  rewritten  to  incorporate  the  artificial 
viscosity  terms,  they  become 


(63) 


3  q  _  1  3  (p  +  Q  +  S) 

cTT  ~  ~  ~PQ  3  x 


Q  +  S) 
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(63) 


v  =  ±_  i2i 

a  Sx 

p  =  f(^  .  U)  =  g(V,U), 

where  we  have  set  the  specific  volume  p  *=V  and  introduced  the  notations 


if  dq/dx  <  0 


if  dq/dx  £  0  . 

Here  i  is  a  constant  having  the  dimensions  of  a  length. 


4  p  dV 

s  =  -  7—  — 
J  v  at 


(64) 


(b)  Difference  Relations 


Let  Ax  and  At  be  small  increments  of  the  variables  x  and  t.  Then  the  set 
of  points  in  x,  t-space  given  by  x=jAx,  t=nAt,  where  j,n=0,  1,2,  •  •  •  is  called 
a  net  (or  grid  or  lattice)  whose  mesh  size  is  determined  by  Ax  and  At.  The 
approximation  to  g(x,  t)  at  the  point  (jAx,  nA  t)  is  denoted  by  g(jAx,  nAt)  or 
simply  by 

(65)  gj°  =  g(j Ax,  nAt)  , 

There  will  be  a  need  to  define  some  of  the  dependent  variables  at  space  mesh 
stations  midway  between  those  corresponding  to  integral  values  of  j.  The 
values  at  these  substations  are  denoted  by 


<66>  «j"t,/2  =  T<Sjn+8j\ 


For  space  differences  the  following  notation  will  be  used 
(6g)jn  =  g((j+  l/2)Ax,  nAt)  -  g  ((j  -  l/2)Ax,  nAt) 

(67)  _  n  n 

“  gj  +  1/2  '  gj  -1/2  . 

Many  finite  difference  systems  have  been  devised  for  problems  involving 
perfect  fluids,  Reference  35.  The  one  most  used  will  be  adjusted  to  accom¬ 
modate  formulation  (63)  and  (64): 
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(68) 


X.n  +  1  .  X  n 


At 


-i-=qn+  1 

J 


(69) 


9jn  +  1  -  q”  1  (6  P)  “  +  (6  Q)jn+  (6  S)  " 


At 


Ax 


(70) 


+  1  = 
j  +  1/2  p 


+  1  „n  4-  1 

1  j  +  1  '  i 


(71)  Un++  *  -  Un  . 

J  +  1/2  j  *  1/2 


At 


A  x 
n 

Pi  ,  i  i2  n  +  1  n  +  1 

J  +  Qj  +  1/2  +  Sj  +  1/2 


+  1 


IVn+  1  vn 

J?  J  J  +  W2  -Vj  +  1/2 

“  f  3  o| 


At 


Vn  +  *  -  Vn 

j  +  1/2  j  +  1/2 

At 


,_0.  n  +  1 

(?2)  P.  =  f 
J  +  1/2 


.n  + 


/n+  1 
j  +  1/2 


-  .  .  11  +  1  TT11  +  1 

Uj  +  1/2  =  g|Vj  +  1/2  '  j  +  1/2 


,73)  ^1/2=- 


8 


Mo 


,n  vn  -  J 

rj  *  l/2  "  Vj  +  1/2 
At 


3  ..n  i  “  1 

Vj  +  i/2  +  Vj  +  1/2 


(74)  Qn  ,  1 

j+1/2  1 


2(apQAx)^ 

0 


fv“  -  vn_  1 
j+1/2  Vj+l/2l 

At 


if  Vn  .  vn  ”  * 
j+1/2  <  * j+1/2 


Otherwise 


In  (74)  we  have  set  l  =  aAx,  where  a  is  a  dimensionless  constant  approxi¬ 
mately  1.  5  to  2,  0, 

Equations  (68)  through  (74)  are  the  formulas  required  to  carry  out  the 
stepwise  computations.  One  starts  the  calculation  by  obtaining  initial  values 
(i.  e, ,  for  time  t  =  0  or  n  =  0)  of  all  quantities  at  all  stations  (and  substations 
j  +  1/2),  The  values  of  the  quantities  for  the  instant  of  time  corresponding 
to  n  =  1  are  first  obtained  at  all  stations.  These  values  are  then  substituted 


25 


into  the  equations  to  evaluate  the  quantities  at  all  stations  for  the  instant  of 
time  corresponding  to  n=  2,  etc.  The  procedure  used  in  going  from  the 


instant  n  to  the  instant 

n  +  1  is  as  follows: 

a) 

Compute 

n  +  1 

q 

from  (69) 

b) 

Compute 

Xn+  1 

from  (68) 

c) 

Compute 

Vn+  1 

from  (70) 

d) 

Compute 

sn  +  1 

from  (73) 

e) 

Compute 

Qn+  1 

from  (74) 

f) 

Compute 

un+  1 

from  (71) 

g) 

Solve  (72)  for 

n  +  1 

P 

In  all  these  calculations  the  values  at  the  substations  j+  1/2  are  computed 
according  to  (67).  A  flow  chart  describing  the  step-by-step  numerical 
procedure  is  given  in  Fig.  16. 

(c)  Initial  and  Boundary  Conditions 

To  use  the  above  scheme  of  calculation  the  initial  values  (i,  e.  at  time 
t*  0  and  at  all  points  of  the  mesh)  of  all  the  dependent  variables  must  be 
known.  At  the  instant  of  impact  their  values  are  known  everywhere  except  at 
the  interface  between  the  impacting  bodies.  These  initial  boundary  values 
may  be  approximated  by  applying  the  Rankine -Hugoniot  relations  to  the  abrupt 
pressure  profiles  which  emanate  from  the  impact  interface. 

The  relations  are  derived  from  the  conditions  that  mass,  momentum  and 
energy  are  conserved  across  the  shocks.  The  corresponding  equations  are, 
Reference  36, 


(75) 


Po  Uo=  P1  U1 


PoUo  +po=  P1  U1  +  P1 
1/2uo2+  Uo  +  PoVo  =  1/2  ul2  +U1  +P1V1  • 


The  equation  of  state  of  the  material  provides  a  fourth  relationship.  The 
subscripts  "0"  and  "1"  refer  to  material  in  front  of  and  behind  the  shock 
wave  respectively.  Here  u.  is  the  flow  velocity  of  the  material  relative  to 
the  shock  wave. 
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Figure  15a  illustrates  the  impact  situation  in  laboratory  coordinates 
before  impact  and  at  the  instant  of  impact.  In  Fig.  17a  the  nomenclature 
used  to  depict  the  situation  just  after  impact  is  illustrated.  Since  the  two 
bodies  may  be  of  different  materials  an  additional  subscript  has  been  intro¬ 
duced  to  make  this  distinction.  The  initial  data  is  represented  schematically 
in  Fig.  17b.  In  Appendix  A  there  are  presented  the  detailed  algebraic  manip¬ 
ulations  necessary  to  obtain  the  initial  boundary  data  when  the  impacting 
bodies  are  of  different  materials. 


In  the  special  case  that  the  two  impacting  bodies  are  of  the  same  material 

d  =  d,  =  d  and  (A-17),  (A- 18)  reduce  to  the  single  equation 
1  ^ 


(76) 


f 


fd 


Now  (76)  can  easily  be  solved  for  d  by  a  machine  program,  and  the  dependent 
variables  subsequently  evaluated  by  (A-9)  through  (A-16), 


pl,(3o  =U  -Po/drl  qi  =  Vo/Z_l 

Wj=  (vq/2)  (2  -  Pj/  P  Q)  (Pt/  PD-  1) 

W2=(vo/2)  (P  i  /  P  D)  (P1/P0- 1)_1 

The  calculations  have  been  carried  out  for  iron-iron,  copper -copper , 
aluminum -aluminum,  cadmium -cadmium,  tin-tin  and  lead -lead  impact.  The 
results  for  p.  vs  v  ,  U  vs  vq  and  P^lp  Q  vs  VQ  are  plotted  in  Figures  18,  19 

and  20  respectively.  In  Figure  21  the  corresponding  Hugoniot  curves  pi  vs 

V  ^  (  =  1  /  p  )  are  plotted. 

(d)  Stability  and  Convergence 


(77) 


p.=  v  d/4 
l  o 

U=  v  2/8 
o 


Having  formulated  the  governing  differential  equations  as  well  as  the 
initial  and  boundary  conditions  in  finite  difference  form,  consider  next  the 
choice  of  space-mesh  size  Ax  and  time  mesh  size  At.  Suppose  that  A  x  has 
been  assigned  a  value  sufficiently  small  to  allow  a  satisfactory  definition  of 
the  flow.  If  At  is  chosen  too  large  the  calculation  will  not  converge  but  will 
oscillate  with  increasing  amplitude  as  n  increases.  This  phenomenon,  called 
instability,  has  nothing  to  do  with  round-off  error,  but  is  a  property  of  the 
difference  system  (68)  through  (74).  In  fact,  the  error  would  only  be  made 
worse  by  taking  a  smaller  value  of  A',  unless  At  is  also  suitably  reduced. 
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In  Appendix  B  a  stability  analysis  of  the  explicit  difference  system  is 
presented  in  some  detail;  the  essential  results  are  contained  in  (B-22)  and 
(B-24). 

To  apply  criterion  (B-22)  we  recall  that 

dV _ 1_  8q_ 

3T  p  9  x 

and  approximate  dq/dx  by-vQ/2^x.  Here  vq/2  represents  the  velocity  of  the 

interface  as  calculated  from  the  Rankine-Hugoniot  conditions.  In  our  calcu¬ 
lations  we  have  chosen  l- a  Ax,  a=  1.5,  so  that  the  criterion  reduces  to 


(78)  **  -<  ,  Po  V 

(Ax)  -  »  +4  5  p  v  Ax 
3  *  o  o 

The  value  for  V=  1/pmay  be  approximated  by  the  Hugoniot  curve  for  the 
particular  material,  Fig.  21. 

Now,  in  regions  away  from  the  shock  p  »  Q  T  and  condition  (B-24) 

becomes  roughly  ° 


At 

Ax 


£ 


/o 

2  V2  I  i^  p !% 


The  adiabatic  sound  speed  c  is  equal  to  * 


c  VP  V  • 
o  »  o 


Consequently,  in  regions  other  than  those  containing  a  shock, 


(79) 


=  0) 


is  the  approximate  stability  criterion. 


*  In  an  adiabatic  process  dU  +  p  dV  =  0,  Since  p  =  g  (V,  U)  we  may  write 
dp  =  dVdg/3V  +  dU  dg/SU.  Elimination  of  dU  between  the  two  yields 
dp  =  (dg/dV  -pdg  /dU)  dV  =  V2  (pdg/dU  -  3g/5V)dp,  whence 

c2=  dp/d/0=V2(p3g/oU-Sg/aV). 
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If  the  same  size  mesh  is  to  be  used  throughout  the  flow  the  more  severe 
restriction  on  At  must  be  observed,  i.  e.  stability  criterion  (78).  Both  sta¬ 
bility  criteria  have  been  calculated  for  each  of  the  parameter  combinations 
considered.  The  results  are  listed  in  Table  I  along  with  the  corresponding 
value  of  A  t  for  the  choice  ^x=  0.  1  cm. 

The  above  conclusions  have  been  substantiated  by  a  number  of  machine 
calculations;  the  results  are  displayed  in  Table  II,  (Reference  37).  The  cal¬ 
culated  requirement  on  At  is  seen  to  be  a  conservative  estimate  in  each  case. 
It  may  l^e  expected  that  the  theoretical  stability  criterion  is  close  to  the  re¬ 
quired  condition  for  all  the  parameter  combinations.  For  the  most  severe 
case,  At  =  5.  1  x  10’  ,  this  means  that  approximately  10,000  cycles  are 
required  for  a  5  microsec.  run.  If  70  space  mesh  points  are  required  this 
means  700,000  point  calculations  which  would  require  approximately  1.4 
hours  on  the  IBM  7090. 
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TABLE  II 


Summary  of  preliminary  runs  testing  stability  of  explicit  difference  scheme 
applied  to  iron-iron  impact.  a  =  1,5  ^  =1.  5  Ax,  Ax  =  0.1. 


NO. 

NET  SIZE 

RESULT 

Theory 

Trial 

1 

0.  100 

Unstable 

1 

0.  03 

0.  050 

Stable 

1 

0.  03 

0.  065 

Stable 

2 

0.  027 

0.  05 

Stable 

3 

0.  014 

0.  02 

Stable 

4 

0.  0023 

0.  005 

Unstable 

4 

0. 0023 

0.  0025 

Unstable 

4 

0.  0023 

0.  002 

Stable 

13 

0.  03 

0.  05 

Stable 

17 

0.  0021 

0.  05 

Unstable 

17 

0. 0021 

0.  025 

Unstable 

17 

0.  0021 

0.  01 

Unstable 

17 

0.  0021 

0.  005 

Stable 
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IMPLICIT  DIFFERENCE  SCHEME 


In  this  section  there  is  outlined  an  alternate  computational  scheme  which 
has  been  developed  to  reduce  the  machine  time  required  for  those  parameter 
combinations  with  large  (1°  and  vq  values.  It  is  an  implicit  scheme  based  on 

Eulerian  formulation,  whereas  the  original  scheme  is  explicit  and  is  based  on 
the  Lagrangian  formulation.  The  new  scheme  is  only  valid  for  impact  between 
bodies  of  identical  material  with  PQ^o.  Therefore  it  does  not  supersede  the 

explicit  scheme,  but  serves  as  a  desirable  complement. 

The  difference  scheme  has  been  adopted  upon  the  suggestion  of  Dr. 

Herbert  Keller,  Institute  of  Mathematical  Sciences,  New  York  University. 

In  treating  similar  systems  of  equations,  Dr.  Keller  has  found  the  implicit 
scheme  to  be  unconditionally  stable.  Thus,  no  restriction  on  At  and  Ax  is 
involved,  only  the  desired  accuracy  need  be  considered  in  choosing  the  incre¬ 
ment  sizes. 

(a)  Difference  Equations 

Since  the  two  impacting  bodies  are  identical  the  phenomena  are  symmetric 
about  the  center  of  mass  coordinates  (Fig.  15b).  The  calculations  will  be 
made  only  for  body  2  where  the  fixed  space  coordinates  are  denoted  by 
(Fig.  22a)  j  =  0,  1,  2,  ...  There  dq/dz£0  and,  consequently,  the  Eulerian  for¬ 
mulations  of  the  governing  equations,  equations  (49)  through  (52),  reduce  to 


(80) 

3 

_£ 

+ 

3(Oq) 

-  n 

dt 

dz 

(81) 

aq 

+ 

q 

aq  _ 

1 

L^  +  ±m0 

at 

dz 

P 

l  dz  3  0 

(82) 

au 

+ 

q 

au  m 

J_ 

M-P+t“ 

3t 

dz 

P 

3<  L  3 

(83) 

p  = 

f 

(P, 

U). 

A  centered  difference  scheme  is  used  to  provide  more  accuracy.  To 
illustrate  the  process  the  details  will  be  carried  out  for  equation  (80).  The 
scheme  is  centered  at  point  z=jAz,  t=(n  +  1/2)  At  a®  shown 
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r  ■ 

h  i  ( 

v  , 

mmi 

i 

O  points  used  in  (80),  (82) 

t 

At 

J _ _ 

1 

1 

1 

) 

y  , 

C 

*  / 

7  •  < 

]  l 
i 

1 

\  1 

"  •  points  used  in  (81) 

X  denotes  center  point 

L.  . 

«-  A  z  -p, 
J 

-  1  j 

i  j- 

hJ 
f  i 

DISTANCE 


in  the  sketch.  Equation  (80)  is  replaced  by 
n+l  n  n  +  1  n 

1  i^l  Pj  -  1  A  1  “Pj  +  1  ( 

2  |  At  St  j 


.  1  J(pCi 1  -(p^).n:  I1  r  H 

+  T  1 - 5“T^ - -  +  - * - =TT- — - >  = 


1  (pq|ntl 

2  |  ‘  '  2Az  | 

Similar  centered  schemes  may  be  written  for  (81)  and  (82).  Upon  simpli¬ 
fication  the  system  becomes  v 


_n  +  1  _  n  +  1  „  n  ,  n 

Pj  -  1  Pj  +  l=Pj  -  1  +Pj  +  1 


(84)  l+U t  j 


(85) 


2  1  [(Pq)  n  +  1  +  (P  q)  n  1-1(0  q)  n  1  +  (pq)  n  1  | 

°  ,L  j  -  1  j  -  1J  L  j  +  1  j+  1-1  i 

/  n+l  n\  r 

(qj  +qj )  b 

f  r_n 


n+1-„nl^/nn+1,n\  L'n  +  l  _  n  +  1  n 


qj  +  8  Az 


*  i -qj  - 1  +  V  i  ■'1 


■j-  i 


=  J.  ^  *  _ \ _  I  I"  n  +  1  n+l  n  n  1 

2  *z  Pjn+  l+P.n  b  j  +  i~pj - 1  +pj  + 1  _pj-i  J 

+  n  +  1  n+1  n  ,  n.  n 

+  3  A  z  Lqj  +  1  2*j  +*j-i  +  V  l_2qj  +  qj-l 
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un+11  +u.n t; 

j  - 1  j  + 1 


-h” . +ujn- .] 

If. [V-l +  V*! [ujn  +  i  "Vm-V- 


n+l  n+l  n  n 

(86)  H  ^ 

(  ’  ~^pn+1  +  pn+1  +  p"  +  p" 


j-1 


j  +  1  j  +  1  j-1 


-2Vrru 


1 

n+l 

n  +  l 

n  n 

2 

n+l  n+l  n  n 

2 

pj  -i 

+  pj  +i 

+  p.  ,  +  p 

>  - 1  j  -  j 

+  3 

Az 

qj  +1  ”  qj  -  1  +qj  +1  "qj  - 1 

The  required  pressure  values  are  computed  from  (83)  according  to 


(87) 


n+l  ,  ,  n  +  1 

Pj  .  1  -M  A. 


un  +  £  \ 

5  -  i  •  Vi  1 


n+l  ,  .  n  +  l  TT  n  +  1  . 
Pj  +i  =  f(Pj  +1  .  u.  +1  ) 


Vl=  £^jn+l*  Ujn+  1> 


n 


Pj  -  i=f(Pj  -!•  Uj  -1>* 

(b)  Initial  and  Symmetry  Conditions 


The  initial,  n  =0,  values  of  the  dependent  variables  p,  q,  U,  p  are  all 
known  from  the  Rankine -Hugoniot  equations  (Fig,  22b).  To  determine  their 
values  at  all  other  time  intervals  a  method  of  computing  the  values  at  time 
n+l  from  those  known  at  time  n  must  be  made  available.  An  iteration 
technique  will  be  devised  for  this  purpose.  The  following  relations,  which 
follow  from  the  symmetry  of  the  problem,  will  be  utilized: 


(88) 


n  n 


n 

q 

-j 


q  n  (hence  q  n  =  0) 

'  j  0 
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(c)  Iteration  Procedure 

To  start  the  iteration  let 

ion\  ~n  +  1  n  ~n  +  1  TTn 
<891  q  i  (o>  =q  j  •  uj(0)=uj 

be  zero-order  approximations  to  q.n  +  ^  and  uj1  +  *  respectively.  Then 

J  * 

the  first-order  approximations  are  calculated  in  the  following  sequence:* 

(a-1)  Substitute  q°.+  1  =  q°.  into  (84)  and  calculate  pn  ^  from  the 
resulting  two-term  recurrence  relation  (n+1  fixed,  j  varied): 

~  n  +  1  ~  n  +  1  n  +  n  1  At  ^  fon ,  t  *  q  n  +  1  +  (oaln  1 

Pj-KD  +  °j  +1(1)  =  P}  -i  Pj  +i  +  7  Az) lj  -MU  qj -KO)  (pq)j  -  lj 


mn\  ~n  +  1  — n  +1  ,  ,n  1 

<90)  qj  ♦  KO)  *  <M'j  +lj 


Several  methods  of  performing  the  calculations  are  discussed  below  in 
section  (d). 

(b-1)  To  obtain  a  trial  value  for  pn.  +  1  use  (87)  with  /txn  +  *  =  P?,t  1  an<l 
TTn+l_~n  +  l  J  J  J'D 

j  j(0)  : 

/on  ~n  +  1  _  f/~n  +1  ~  n  +  1\ 

(91)  p  j(  lM  .  Uj(0)  j. 

(c-1)  To  obtain  a  first-order  approximation  for  U™  +  *  use  (86)  with 

9j  .  Pj  .  Pj  replaced  by  qj  ,  Pj(J)  »  Pj  {1.)  respectively: 

ft  n+1  ~  n+1  T.  n  -] 

j  -1  (1)  J  +1  (1)  L  J-l  Uj+l 

[V-«'>  *  Vuo,%;,  *  'j-ijfSjr.lt, 

~>n  +1  —  n+1  n  n  . 

(92)  _£t  qj  +  l  (0)  ~qj -1(0)  +  qj  +  l  ~qj-l  f  _  . 

“  A~z  Z~HTl - ~~ZrT+  i - ~n - )  -  2VT  To  (continued) 

Pj-  KD  Pj  +  1  (1)  +  Pj+1  +Pj-1  * 


(continued) 


*  This  procedure  is  adopted  upon  the  suggestion  of  Dr.  Herbert  Keller, 
Institute  of  MathematicalSciences,  New  York  University. 


34 


1  L  n  +  1  .  -  n  +1  n,  n  1 

"I  [P  j-  l(l/)+  Pj+1  (1')  +  Pj+l  +  Pj-lJ 

2  u  f  ~n+l  -n  + 1  n  n  "I  | 

T  13 +i  (0)  "  j_1(0)  qj+1  ‘  qi  -IJ)  • 


As  with  p,  we  have  here  a  two-term  linear  recurrence  relation  for  . 

Its  solution  is  also  discussed  in  section  (d).  ' 

(d-1)  A  first-order  approximation  for  p^11"1"*  i®  now  computed  from  (87)  with 
n  +1  -n  +  1  and  T7n+1  ~n  +1 

pj  -°j(ll  Vi  '“jUI 


/r\o\  ^n+1  cl +*  n+1  Cln+1  \ 

(”1  pi(D*  Tim-  uj(i)j- 


n+1 


(e-1)  To  calculate  a  first-order  approximation  for  use  (85)  with 

n+1  —  n+1  n+1  —  n+1  ,  ,  n+1  .  ... 

p j  =  p.  ,  p^  =  Pj  (j)  and  replace  in  the  difference  expression 


for  q  3  q/  3  z  by  q  . 


-n  +  1 

q j  ( !)  * giving 


j  (0) 


The  latter  substitution  linearizes  the  equation  for 


(94) 


.  1  At 


-n+1  n  1  At  /-n+1  ,  n\  f— n+1  ~  n +1  ,  n  n  1 

qj(l)-qj+8  at  (qj(0)+  qj  ;[qJ+M0)  -qj  -  K0)+qj+l  -  qj  -  ij 

1  l^n+l  —n+1  _n  _  n  "| 

n+1  n)  pj  +1(1)  ”Pj  -  1(1)  pj  +  1-  '  Pj  -  1 

j  (1)  P  j  lL  J 

n  o  n  n  V 

■j  +i  -2,j  *  ’j-ijj. 


2  Az  ~ 


*  3  A  z 


-n+1  n  +  1  .-n+1 

qj  *iui  -2<1j  id  *  ^  - 


This  is  a  three-term  linear  recurrence  relation  for  q”  jjj  (n  +  1  fixed,  j 
varied);  a  method  of  solution  is  discussed  below  in  section  (e). 


The  second-order  approximations  are  calculated  by  merely  repeating  the 
iteration  process:  (a-2),  (b-2),  ...»  (e-2).  The  resulting  equations  differ 
from  the  corresponding  equations  of  the  first  iteration  only  in  that  the  sub¬ 
scripts  (1),  (2)  replace  the  subscripts  (0),  (1)  respectively. 

In  general,  to  proceed  from  the  k-order  approximation  to  the  (k+1)  -  order 
approximation  we  go  through  the  above  iteration  process  (with  subscripts  (0), 
(1)  replaced  by  k  and  k+1  respectively).  The  process  is  repeated  until  a  rea¬ 
sonable  convergence  criterion  is  satisfied.  Usually,  only  a  few  cycles  are 
required  in  such  schemes. 
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Assume  that  it  has  been  decided  that  K  iterations  are  sufficient.  Then 


set 


„n  +  1  ~  n  *•  1 

p.  •=  p. 


n  +•  1  —  n  +  1 

J  (K)  *  qj  =  qj 


TT  n  f  1  rfnfl  nU  ^n*  1 
(K)  '  Uj  "  j  (K)  *  pj  (K)=  Pj  (K)* 


and  proceed  to  the  n  +  2  time  step.  A  flow  chart  describing  the  numerical 
procedure  is  given  in  Fig.  23. 

(d)  Two-Term  Recurrence  Relations 

It  is  a  consequence  of  the  symmetry  relations  (88)  that  when  j  =  0  the  two- 
term  recurrence  relation  (a-1)  simplifies  and  yields  the  explicit  formula 


I  e  \  n  +  1  n 

(95)  #m)  =  Pi 


.  1  At  n 

1  ~  2  Az  q! 

,  .  J_  At  ♦  1 
1  *  2  Az  ql  (0) 


Now  (a-1)  may  be  used  to  compute  p\n  ^  at  all  the  odd  j  space-mesh  points: 


(96) 


n<-  1  f.  1  At  ^n+’ll  ^ntlf.  1  At  „n  ♦  H 

Pj  +2(l)[_1+2  Az  qj^2(°)J+Pj(l)  ‘  2  Az  qjd)  J 

nn  Ti  1  Al  n  ]  an  f,  2_At  nl 

“°jt2  L1  -2  4.^  +  2]+^  [l*  2  Az  1j  J 

Thus,  set  j  —  1  and  compute  p^n  *  ^  in  terms  of  Pj11  ,  set  j  —  3  and 

,  ~  n  ♦  1  .  .  ,~n+l 

compute  p,.  (l)  ln  terms  of  P3  (i)  *  ®tc* 

n  t  \ 

The  /ST .  ^  at  even  j  are  determined  by  the  continuity  of  p - for 

n  4  1' 

sufficiently  large  j,p^  =  pQ  (the  density  of  the  undisturbed  medium, 
which  has  not  yet  been  reached  by  the  shock  wave).  Let  J  denote  such  a 
large  even  integer  at  time  t  =  (n  ♦  1)  At.  Then  J5j  ~  PQ  an<*  t^ie  recurs*on 
formula  is  used  to  calculate  from  right  to  left.  Thus,  set  j  =  J  -  2  and 
compute  p  j  2  ( 1)  in  terms  °f  P  j  ^«setj  =  J-  4  and  compute  p  j  4  ( i) 

,  r-n  +  1 

in  terms  of  pj  2(1)'  etc* 


A  suitable  value  for  J  may  be  found  by  first  making  the  calculations  for 
the  odd  valued  mesh  points  until  p*1  has  decreased  to  the  value  pQ .  This 


value  is  then  taken  to  be  J  -  1.  This  left-to-right-to-left  technique  is  depicted 
in  the  sketch. 
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NUMERICAL  RESULTS 


The  various  combinations  of  assumptions  for  Tq  and  vq  for  which 

calculations  have  been  made  are  displayed  in  Table  I.  In  computing  the  values 

for  dimensionless  parameters,  B  and  R  ,  shown  there  the  characteristic 

o  o 

length  is  taken  as  unity.  This  choice  is  meaningful  only  after  a  stable  wave 
has  been  established.  In  Figures  24  through  35  the  results  of  the  calculations 
are  illustrated  by  pressure  profiles  for  a  number  of  typical  parameter  com¬ 
binations.  Some  additional  profile  plots  have  been  presented  elsewhere 
(Reference  38). 

To  check  the  accuracy  of  the  program  the  calculations  for  cases 

No.  1  (v  =  0.  5,  T  =  0,  a  =0)  and  No.  49  (v  =  1.0,  T  =  0,  u  =  0)  have 
o  o  o  o  o  o 

been  examined  in  detail.  The  computed  pressure  profiles  are  depicted  in 

Figures  24  and  34,  The  figures  show  that  a  stable  shock  front  is  quickly 

established.  For  comparison,  the  Rankine-Hugoniot  solutions,  applicable 

since  here  u  =  r  =0,  are  also  shown  at  t  =  1  and  t=  5  microseconds,  The 
Ko  o 

other  calculated  dependent  variables  behind  the  stable  shock  may  also  be 
compared  with  the  corresponding  Rankine-Hugoniot  values  for  those  cases 
where  U0=T0=  0.  This  is  done  in  Table  III  for  cases  No.  17  (vQ=  4.0,  T0=  0,  u0=  0)  and 

No.  33  (vq  =  7.  5,  Tq  =  0,  Mo  =  0)  as  well  as  Nos.  1  and  49.  The  computed  quantities 

represent  mean  values  about  which  there  are  small  oscillations  at  the  various  mesh 
points  behind  the  front.  The  agreement  is  seen  to  be  quite  satisfactory. 

The  effect  of  viscosity  on  the  amplitude  and  duration  of  the  transient 
pressure  profile  which  occurs  immediately  after  impact  is  illustrated  by 
comparing  Figures  24a,  25a,  26a  and  27.  For  all  these  cases.  Nos.  1,  2, 

3  and  4  respectively,  the  impact  velocity  and  yield  strength  are  the  same 
(vq=  0.  5,  Tq  =  0).  When  M  Q-  0  (No.  1)  the  interfacial  pressure  builds  up 

to  the  final  stable  profile  amplitude.  For  fx  q>0,  however,  the  material 

near  the  interface  is  seen  to  be  subjected  to  very  large  pressures  in  the  first 

few  tenths  of  microseconds.  Subsequently,  the  interfacial  pressure  decreases 

from  the  original  higher  magnitude  to  the  final  stable  profile  amplitude.  For 

u  =  0.  08  (No.  2)  the  duration  of  the  overshoot  is  less  than  0.4  microsecond; 
r  o 

for  u  =  0.  8  (No.  3)  the  duration  is  less  than  0.8  microsecond;  for  a  =8.0  the 
o  o 

duration  of  the  overshoot  persists  for  5  microseconds.  The  amplitude  of  the 

transient  pressure  pulse  also  increases  with  an  increase  in/i  . 

The  effect  of  viscosity  on  the  characteristics  of  the  stable  pressure 
profile  is  illustrated  by  comparing  Figures  24b,  25b,  26b  and  27.  These 
again  correspond  to  cases  Nos.  1,2,3  and  4  respectively.  When  /i  =  0 

andMo  :0.08  the  thickness  of  the  stable  profile  is  about  0.2  cm.  This  is  due 
to  the  artificial  viscosity  term.  But  for  the  thickness  becomes  0.  7  cm 
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TABLE  III 


Rankine -Hugoniot  values  compared  with  finite -difference  calculations  for  one¬ 
dimensional  iron-iron  impact  (Mc“*  ^  “  0)*  VP  denotes  velocity  of  stable 
profile  relative  to  the  interface.  0 


Case  1  | 

Hugoniot 

P 

1.  565 

1/P 

0. 0873 

U 

0.  0312 

VP 

0.  5466 

<vo-  •  5>  < 

[  Computed 

1.  541 

0.  0867 

0. 0265 

0.  534 

Case  49  | 

|  Hugoniot 

4.  646 

0.  0734 

0.  125 

0.  682 

<vo=  1 

[  Computed 

4.  601 

0. 0734 

0.  122 

0.  673 

Case  17  I 

^  Hugoniot 

49.62 

. 04662 

2.0 

1.  157 

II 

0 

> 

1  Computed 

48. 5927 

.  046869 

2. 02275 

1.  159 

Case  33 

^  Hugoniot 

154.2 

.  03601 

7.  031 

1.481 

(vQ=  7.5)1 

[  Computed 

153. 6917 

.  036368 

7. 1429 

1.  462 
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and  when  a.  =  8.  0  it  is  greater  than  5  cm.  In  the  latter  case.  R  &  B  <1 
'O  o  o 

and  it  might  be  anticipated  that  the  viscous  terms  would  predominate.  The 

amplitude  and  velocity  are  only  slightly  increased  by  increasing  /i0.  The 

time  required  for  the  stable  profile  to  be  established  varies  from  less  than 

2  microseconds  for  q=  0  to  approximately  10  microseconds  for  8.0, 

The  effect  of  impact  velocity  is  illustrated  by  comparing  Figures  27,  35 
and  31.  These  correspond  to  cases  No.  4  (v  =  0.  5,  Tq=  0,  jzq=8),  No.  52 

(v  =1.0,  T  =0,  n  =8)  and  No.  20  (v  =  4.  0,  T  =  0,  a  =  8}  respectively.  The 
O  Of.  o  o  o  o 

duration  of  the  transient  pressure  pulse  is  observed  to  decrease  from  about 

5  microseconds  for  v  =0.  5  to  about  1  microsecond  at  v  =1.  0,  to  less  than 

o  o 

0.6  microsecond  at  vq=4,  0.  The  stable  profile  becomes  steeper  as  the 

impact  velocity  is  increased.  Naturally  the  amplitude  and  velocity  of  the 
wave  are  increased  when  v  is  increased. 


To  illustrate  the  effects  of  varying  nQ  and  ,  results  for  cases  where 
To=0  have  been  cited,  i.  e.  results  for  viscous  liquids;  the  same  effects  are 

produced  if  TQ  is  held  equal  to  some  reasonable  positive  value.  To  under¬ 
stand  the  behavior  of  such  a  visco-plastic  medium  refer  to  equation  (62)  and 
observe  that  in  front  of  the  disturbance  dq/dx=0,  in  the  disturbance  dq/dx<0 
and,  finally  3q/}x  drops  back  to  zero  after  the  disturbance  passes.  Conse- 


2  2  2 
quently,  r  >  Tq  in  the  disturbance,  but  T  drops  to 

means,  in  terms  of  our  visco-plastic  model,  that  there  is  flow  only  in  that 
part  of  the  medium  through  which  the  disturbance  is  currently  passing;  it 
again  becomes  rigid  behind  the  disturbance.  Thisphenomenonis  illustrated 


Tq  after  it  passes. 


This 


2 

in  Figure  36  where  the  value  of  T  (which  is  a  measure  of  the  distortion  of 
the  medium)  is  shown  at  various  time  intervals  for  two  typical  parameter 


?  2  . 

combinations  (Nos.  14  and  15).  At  each  instant  T  >TQ  is  seen  only  in  a  finite 
region  which  represents  the  current  position  of  the  disturbance.  Only  in  this 
moving  region  of  disturbance  does  the  medium  behave  as  a  viscous  liquid. 

The  region  is  more  spread  out  the  greater  the  value  of 


Comparison  of  Figures  26b,  28  and  29  (casesNos.  3,  7,  11  respectively) 
shows  that  even  with  an  impact  velocity  as  low  as  vq=  0.  5  cm/microsecond 

there  is  little  change  effected  on  the  stable  profile  by  inclusion  of  the  yield 

stress  if  it  is  as  low  as  t  =  0.  01  or  T  =0.  1.  Only  small  increases  in  the 

o  o 

pressure  and  in  the  disturbance  velocity  are  apparent.  This  might  be  ex¬ 
pected  since  for  all  these  parameter  combinations  the  ratios  B^:  Rq  and  1;Rq 

are  small,  i.  e.  ,  the  inertial  terms  predominate.  A  condition  to  keep  in 
mind  is  that  in  an  actual  cratering  process  the  flow  velocity  must  always 
decrease  to  the  point  where  t  is  important. 
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On  the  other  hand,  comparison  of  Figure  26  (case  No.  3)  with  Figure 
30  (case  No.  15)  shows  that  both  the  amplitude  and  the  velocity  of  the  pressure 
pulse  are  significantly  increased  by  the  inclusion  of  the  strength  term  if  it  is 
as  large  as  tq=  1.0.  The  shape  of  the  pulse,  however,  is  apparently  not 

strongly  affected,  nor  is  the  tim^  required  for  the  profile  to  be  stabilized. 

For  these  cases  the  inertial  and  trength  terms  are  both  important, 

R  :  B  =  1. 97. 
o  o 

The  characteristics  of  thi  stable  profiles  for  all  the  various  cases  are 
depicted  in  Figures  37,  38  and  39.  There  p  denotes  the  pressure  behind  the 
disturbance;  AP  and  VP  denote  the  thickness  and  velocity  (relative  to  the 
interface)  of  the  stable  pressure  profile,  respectively. 
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REMARKS  ON  THE  EQUATION  OF  STATE 


In  our  model  the  medium  is  considered  to  act  as  a  viscous  liquid  rather 

2  2 

than  a  solid,  provided  r  >Tq  .  It  is  therefore  consistent  with  ordinary  viscous 

fluid  theory  that  the  equation  of  state  be  given  in  terms  of  a  hydraulic  pressure 
p,  equal  for  all  directions.  We  have,  as  usual,  taken  p  =  -  (T  +  TgQ  +  rzz)/3 
as  the  thermodynamic  pressure. 


The  particular  equation  of  state  employed  in  the  calculations,  p  =  f  (p,  U), 

was  determined  by  the  Los  Alamos  group  from  measurements  on  pressure 

pulses  induced  by  high  explosive.  The  method  is  indirect  in  that  the  observed 

quantities  are  the  pulse  velocity  and  the  free-surface  velocity  produced  by 

normal  reflection  of  the  pulse  from  a  free  boundary.  Pressure  (strictly, 

stress  normal  to  the  wave  front,  -  r  )  and  corresponding  values  of  internal 

z  z 

energy  and  density  were  computed  from  these  measurements  by  means  of  the 
Rankine-Hugoniot  relations.  An  equation  of  state  based  on  the  assumption  that 
HQ  —  Tq=  0  has  thus  been  employed  to  calculate  the  behavior  of  a  model  for  the 

material  which  assumes  that  these  parameters  are  not  zero.  This  certainly 
leads  to  errors  but  they  are  of  second  order  and  would  not  be  expected  to 
mask  the  effect  of  including  the  viscous  and  plastic  terms  in  the  other  equa¬ 
tions  governing  the  model.  The  results  have  borne  this  out  since,  as  physical 
reasoning  would  imply,  p  chiefly  affects  the  shape  of  the  stable  disturbance 

and  T  has  its  main  effect  on  its  amplitude.  At  lower  velocities,  t  would 
o  o 

also  affect  the  shape. 

Other  remarks  on  the  equation  of  state  are  also  relevant.  In  converting 
the  measured  velocities  to  pressure-energy-density  states  it  was  tacitly 
assumed  that  a  stable,  abrupt  disturbance  was  obtained.  Verification  by 
direct  measurement  of  the  pressure  profile  has  not  been  possible,  and  justi¬ 
fication  for  the  assumptions  is  based  on  the  reasonable  agreement  with  extrap¬ 
olation  of  hydrostatic  data.  This  should  not  be  construed  as  proof  that  the 
viscosity  and  strength  effects  are  negligible,  however,  since  even  with  vis¬ 
cosity  factors  as  large  as  jiQ=  0,8  and  yield  stress  as  great  as  to=0.  1,  when 
Vq  =  .5,  the differencesin  the  velocity  and  amplitude  of  the  stable  pressure 

would  be  difficult  to  observe  by  such  measurements. 

The  fact  that  -  T  and  not  p  is  the  actual  pressure  reflected  from  the 
zz 

free  surface  in  the  experiments  may  not  greatly  affect  the  equation  of  state 
calculations  since  the  shape,  velocity,  and  amplitude  of  the  stable  p  and  -r 

ZZ 

profiles  are  nearly  identical.  To  see  this  compare  Figures  26b,  27  and  30b 
with  Figures  40,  41  and  42  respectively;  these  show  the  two  corresponding 
profiles  for  three  typical  cases.  Nos.  3,4,  15. 

Prior  to  the  establishment  of  a  stable  profile,  however,  the  components 

of  the  deformation  stress  tensor,  of  which  T  +  p  is  one,  are  not  small.  At 

z  z 
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the  instant  of  impact  between  the  two  bodies  a  very  large  velocity  gradient  is 
imposed  on  the  material  near  the  interface.  As  the  front  of  the  disturbance 
propagates  into  the  body  the  gradient  at  the  interface  decreases,  and  the 
gradient  at  the  front  of  the  disturbance  also  decreases  because  of  the  smear¬ 
ing  action  of  the  viscosity.  From  (48)  is  seen  that  the  von  Mises  statistic, 

2 

r  ,  which  is  a  measure  of  the  magnitude  of  the  components  of  the  deformation 
stress  tensor,  must  act  similarly.  Thus  when  viscosity  is  present  the 
material  near  the  interface  is  subjected  to  a  much  greater  distortion  than 
material  away  from  the  interface.  This  is  illustrated  in  Figure  36  by  the 
2 

envelopes  of  the  r  distributions. 

These  latter  observations  are  consistent  with  experimental  evidence  that 
internal  structural  changes  in  a  metal  can  be  related  to  the  distribution  of 
stress  that  existed  in  an  impulsively  loaded  body  by  plotting  contours  of  equal 
hardness  on  sections  of  the  body,  Ref.  39.  Contours  were  found  to  coincide 
with  the  isochromatics  obtained  in  photoelastic  studies,  i.  e. ,  the  contours 
lie  along  lines  of  maximum  shear  stress.  Recent  microhardness  studies 
of  one  dimensional  impacting  plates  have  shown  that  indeed  the  microhardness 
near  the  impact  interface  is  maximum,  the  value  decreasing  rapidly  outside 
the  interfacial  zone.  Ref.  40. 
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CONCLUSIONS 


A  visco-plastic  model  for  hypervelocity  impact  has  been  proposed  which 
takes  into  account  the  inertial,  viscous,  and  plastic  effects.  This  was  ac¬ 
complished  by  introducing  a  viscosity  factor  M0  and  a  dynamic  yield  stress  TQ 
into  the  perfect  fluid  equations.  From  an  examination  of  the  resulting  system 
of  equations  several  dimensionless  parameters  were  found  which  control  the 
relative  importance  of  the  three  effects  at  the  various  stages  of  the  cratering 
process.  The  inertial  effect  was  found  to  be  important  throughout  the  early 
stages  while  the  strength  of  the  medium  is  dominant  during  the  final  stages. 
Immediately  after  impact  the  viscous  effect  is  very  large  in  the  zone  near  the 
contact  interface.  Its  magnitude  decreases  as  the  strain-rate  gradient  decreases, 
but  it  may  remain  important  throughout  the  flow  process.  The  viscosity  also 
has  the  important  function  of  introducing  anisotropic  stresses  into  the  flow 
which  are  essential  if  the  strength  effect  is  to  be  included. 

In  the  absence  of  definitive  data  for  MQ  and  T0  in  the  hypervelocity  impact 
regime,  computations  were  performed  on  a  one -dimensional  model  in  which 
the  values  of  these  two  parameters  were  varied.  The  above  qualitative  con¬ 
clusions  were  verified.  Specifically,  the  following  was  found. 

1)  The  assumption  UQ  >  0  results  in  large  initial  values  for  the  pressure 
and  deformation  in  a  zone  near  the  impact  interface.  As  UQ  is  increased  the 
effect  becomes  greater  and  the  disturbance  propagates  a  greater  distance  be¬ 
fore  reducing  to  its  stabilized  shape  and  amplitude.  For  impact  velocity 

vQ  =  0.  5  cm /Msec.  ,  the  time  the  disturbance  propagates  before  a  stable  pro¬ 
file  is  obtained  varies  from  about  2  to  3  microseconds  for  M0  =  .  08  to  more 
than  10  microseconds  for  Mo  ~  8.  The  required  time  is  less  the  greater  the 
impact  velocity;  the  value  of  Tq  has  little  effect. 

2)  The  amplitude  and  velocity  of  the  stable  pressure  profile  are  only 
slightly  increased  as  MQ  is  increased,  but  its  width  (shape )  is  significantly 
larger.  Increasing  the  yield  strength  has  little  effect  on  the  shape  of  the 
stable  pressure  wave;  it  significantly  increases  its  amplitude  and  velocity 
only  if  Tq  is  as  large  as  one  megabar.  The  latter  conclusion  is  valid  for 
particle  velocities  of  0.25,  0.  5  cm /microsecond  and  larger.  At  lower 
velocities  tq  has  a  more  significant  effect  on  the  viscosity  coefficient,  see 
(8),  and  thus  more  effect  on  the  pressure  wave. 

These  conclusions  may  be  related  to  the  qualitative  model  of  crater 
formation  that  has  evolved  from  experimental  studies  in  which  the  actual 
cratering  process  has  been  monitored,  Ref.  41.  The  discovery  was  made  that 
though  only  five  to  ten  microseconds  are  required  to  use  up  the  projectile, 
the  crater  continues  to  enlarge  for  several  hundred  microseconds.  The 
mechanism  of  crater  formation  is  therefore  essentially  one  of  cavitation,  the 
size  and  shape  of  the  final  crater  being  determined  by  (a)  the  shape  and 
amplitude  of  the  pressure  wave  established  during  the  first  five  to  ten  micro¬ 
seconds  by  the  action  of  the  projectile  on  the  target,  and  (b)  the  resistance  of 
the  target  material  to  flow,  (c)  The  flow  continues  until  the  amplitude  of  the 
wave  decreases  below  the  intrinsic  yield  strength  of  the  material. 
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The  calculations  presented  in  this  paper  show  that  the  shape  and  amplitude 
of  the  pressure  wave,  (a),  are  in  turn  strongly  dependent  on  the  viscosity  of  the 
medium.  This  is  especially  true  during  the  first  microsecond  after  impact  when 
the  strain-rate  gradient  is  largest.  Also,  the  resistance  of  the  target  material 
to  flow,  (b),  depends  on  the  viscosity  factor  HQ  and,  to  a  lesser  degree,  on  the 
strength  factor  T Q;  the  viscosity  coefficient  becomes  larger  and  more  dependent 
on  Tq  at  the  smaller  strain-rates.  Finally,  (c),  the  strength  factor  TQ  controls 
the  instant  when  the  flow  ceases. 

Thus,  buth  UQ  and  tq  are  important  in  determining  how  long  the  crater 
continues  to  expand.  This  may  explain  why  a  crater  in  Lucite  stops  expanding 
earlier  than  one  in  aluminum  despite  the  relative  magnitudes  of  their  yield 
strengths. 

To  establish  the  validity  of  this  hypothesis  will  require  experimental 
data  which  are  currently  not  available.  The  necessary  definitive  experiments 
for  the  evaluation  of  UQ  and  T Q  should  be  performed.  They  are  likely  to  be 
one-dimensional  in  character.  The  dependence  of  both  on  the  temperature  of 
the  medium  should  be  studied,  since  this  effect  must  eventually  be  incorporated 
into  the  refined  visco-plastic  model.  Meanwhile,  the  theoretical  program  may 
be  pursued  by  developing  a  numerical  scheme  for  the  calculation  of  the  visco¬ 
plastic  model  under  axial  symmetric  impact  conditions. 
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APPENDIX  A 


In  Fig.  16a  the  notations  used  to  describe  the  wave  fronts  in  the  impacting 
materials  are  depicted.  The  particle  velocities  qm^  and  wave  velocities 

are  with  respect  to  fixed  laboratory  coordinates.  The  pressure  and  particle 
velocity  must  be  continuous  across  the  interface  so  that 


Pll 

-P21  -  Pi 

u10  =vo+Wl 

u20  ~ 

qll 

=  q21  =  qi 

ull=  qi+  W1 

U21  = 

*i  -  "2 

and  application  of  equations  (75)  to  the  two  discontinuity  surfaces  yields 
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To  these  six  equations  in  the  eight  unknowns  p^,  q^.,  Wj,  W^,  Pjj*  ^21'  ^1* 
and  we  may  append  two  more. 


(A-7) 

pi  _fl  (P11* 

Ul> 

(A-8) 

Pi  =  f2  <p2i» 

U2> 

the  equations  of  state  for  the  two  materials. 

Equations  (A-l),  (A-3)  may  be  solved  for  Wj,  W2  and  substituted  into  the 

remaining  six.  The  new  equations  corresponding  to  (A-2),  (A- 5)  may  then  be 
solved  for  p^,  q..  These  manipulations  give  the  following  relations: 

(A-9)  (»„/P10- 1)-1 
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(A- 10)  W2=(q iP21/P2Q)  (P21/P20-1)-1 

(A-ID  Pll/P10  =(1  "  P10/dl)"1  where  dl  S  P11  (P11/P10‘1)"1 


(A-12)  P21/P20  ={1  -p20/d2)_1  where  d2  =  D2l  (P21/P2 ^l)*1 

(A.  13)  qj 

(a- i4)  ^  =v92  d1  d2^1+yi2)-2 . 

Equations  (A-3),  (A-6)  may  be  solved  for  U.,  U«  and  the  above  relations  used 
to  give  1  “ 

(A.  15,  2(fc+fc)-Z 

(A-16)  U2  =i  vo2  '2. 


If  (A-ll),  (A-12),  (A-15)  and  (A-16)  are  substituted  into  (A-7),  (A-8)  there 
finally  result  two  equations  in  the  two  unknowns  dj,  d2: 


(A- 17) 


^10 


(A -18) 


20 


1  -  P20/d2 


Once  d j  and  d2  have  been  determined  the  quantities  q.,  pj,  p^,  /D^,  Uj,  U2 

may  be  computed  from  equations  (A- 11)  through  (A-16).  The  required  initial 
data  are  then  available  for  the  finite -difference  calculations  as  displayed  in 
Fig.  16b. 


Since  the  functions  f  j,  f2  are  too  cumbersome  for  explicit  solution  for  d^ 
or  d2  to  be  practical,  a  numerical  method  must  be  used.  One  such  method  of 
successive  approximation  is  as  follows: 


a)  Guess  a  value  for  d.,  e.  g.  d 


(1) 
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b)  Substitute  ^  into  (A.  17)  and  compute  the  corresponding  estimate  of 

^2*  say  by  trial-and-error. 

c)  Substitute  dj^,  d2^  *n*°  *ke  two  *ides  of  (A- 18). 

(i)  If  equality  holds  dj  =d^*\  d^  =  dg^  and  process  is  completed. 

(ii)  If  equality  does  not  hold  dj  ^d^1*,  d?  t  d2(1)  and  a  second  guess 
for  dj,  say  d^  \  must  be  made  and  the  process  repeated,  etc. 

The  method  is  simple  but  involves  extensive  calculations.  A  program  has 
therefore  been  written  to  carry  out  the  computations  on  a  high  speed  digital 

COmnnfpt*  0  r  © 
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(B-l) 
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) 


where  we  have  eliminated  X  from  (68)  and  (70)  to  obtain  the  second  equation 
and  the  equation  of  state  is  rewritten  as  p  *  f(V"\  U)  =  g(V,  U).  Here  l  is  a 
parameter  with  dimension  of  length  which  essentially  determines  the  magnitude 
of  the  pseudo-viscosity. 


The  analysis  of  the  stability  of  this  system  follows  the  method  outlined  by 
Richtmyer  (Reference  35).  The  equations  of  first  variation  of  (B-l)  will  be 
obtained  in  which  quantities  of  the  second  and  higher  order  are  dropped. 

This  will  give  us  linear  equations  for  the  first  order  variations  4.  V,  S,  0,  £), 
p  (the  dot  does  not  indicate  time  derivatives)  in  which  the  zero  order  quantities 
appear  as  coefficients.  The  equations  obtained  are 
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The  zero  order  quantities  are  considered  constants  and,  consequently, 
superscripts  and  subscripts  denoting  net  points  are  omitted  from  them. 


The  first  order  quantities  are  assumed  to  have  the  Fourier  representations 
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where  x  =  jAx.  The  von  Neuman  stability  criterion  is,  essentially,  that  the 

coefficients  a£ . remain  bounded  as  the  calculations  proceed  from 

t  s  nit  to  t  =  (n+1)  At,  t  =  {n+2)At,  etc.  To  investigate  this  substitute 
representations  (B-3)  into  (B-2)  and  set  like  harmonics  equal  to  zero  to 
obtain  the  relations 
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(B- 10) 


A”k,1  =  Ak-iS4t(FktEktck) 
B”+‘=  b”  4  iSM  a£+‘ 

„n+l 


.n+1 


4  Uo 

f sn+1  -a"  ] 

3  VAt 

L  k  k  J 

•”-(p+QfS+^lTo)  [Bk+1 

-<] 

(P002 

(ft  KH : 

.  l) 

,„2 

J  VAt 

av  rwn+l  n-j 

at  LBk  "Bkj 


Fn+1=l&  Bn+1+±£.  Dn+1 

*k  av  k  a  u  k 


_  2  sin(k  Ax/ 2) 

p  A  x 
o 


Relations  (B-5),  (B-6)  may  be  combined  to  yield 

c”  =  “  4  17  i®Av 

k  3  V  k 

This  may  be  substituted  into  (B-4)  to  obtain 

,B-n>  =  (i  -i  %  ftt)  ^-is«  [e”  +  r”  ] 

If  (B-ll)  is  substituted  into  (B-5)the  result  is 

\K 


(B- 12)  B”+1  =  b£  +  i  0  At  (l  32At)  A”  +  6Z(At)2  [eJ  +  f"  ] 

Now  (B- 12)  may  be  used  to  eliminate  b£+1  from  (B-7)  and  (B-8)  with  the  results 
(B-U)  D^1.d“-(P40.S41(|  tJ  [iSA»(l-i-^-S2At)^  4  92(At)2(E^4F”)] 
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(It)2  K+iS1‘  (>  -1 T  -'K*  sW(e^)] 


*vfr  ttH84'!1  •i?8!4K4!!|4tl\Ek,Fk)  ] 

Finally,  substitution  of  (B-12),  (B-14)  into  (B-9)  yields 

Fkntl'lv{Bk  +  i8"(‘‘f  T02“)AktB2l4«2(Ek+f”)  } 

tfuK'(p+Qts+^F  To)[ie4t(‘-TT  ^‘K^'^K^k)]} 


(B-15) 


Upon  introduction  of  the  notations 


a « c/2)  v  - 1  4  fi2A. 

PqAx  ^  '  “j  y  3  At 


(B-16)  A  3  p+Q+S  +J4  t 

!  3  o 

av  Aau 

e=  -L 

‘»o«*  8V 

2v  at 

w  V  at 

the  relations  (B-H)  through  (B-15)  may  be  written  in  the  matrix  form 
(B-17)  Yn+1  =  GYn 


where 


i 

Cc* 

i _ 

An+1 

Ak 

* 

Bk*‘ 

Dk 

Yn+1> 

°r 

■E 

n+l 

Ek 

Fk 

-,n+l 

Fk 
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and  the  "amplification  matrix",  which  depends  on  and  t  ,  is  given  by 


(B-18) 


0  -10  At 


igyAt  1  0  ^(At)2  p2(At)2 

-ipVAAt  0  1  -^A(At)2  -p2A(At)2 

iSy£4t(|r-c)  -2eC  0  i?CW2(|--c)  -«) 


2  2 
-0A(At) 


ipyAAt 


la  la  s2 


av  au 


P  6(At) 


^6(At>2 


Now,  expansion  of  the  determinant  | G  —  X 1 1  and  setting  the  result  equal 
to  zero  shows  that  the  eigenvalues,  X,  of  G  satisfy  the  equation 

(B-19)  X2(X-1)  |(X-1)2-(X-1)  ^At[fiAt  +  C(2-eAt)-|  -^  ]  +  p2(At)2(2eC- 6)}  *  0 

2  2 

If  Is  constant  and  At/{Ax)  =  0(1)  as  At,  Ax—^O,  then  P  At  a  0(1)  and  the 
secular  equation  reduces  to 

(B-20)  X2(X-l)2|x-l-p2At[2C-  |  -^]-0(At)} 

Thus,  the  von  Neuman  requirement  for  stability 

(B-21)  j  Xm|  s  1+0  (At) 


is  satisfied  provided 


r  2  _^o_  Vi  .av. 
L  3  V  V  I  at 


fl] 


where  we  have  used  the  fact  that  3V/at  <  0.  The  inequality  will  hold  provided 


(B-22) 


«_  %y 
..  ,2  6  . 


(Ax)'1  f  M0+4(po4)2||^| 


It  may  be  noted  from  (B-22)  that  in  the  limit  as  the  viscosity  tends  to 
zero  the  stability  criterion  reduces  to 


(B-23) 


V 

2  i  a  V 


(Ax)  4^|ff 


K  ■  0) 


which  is  the  value  given  by  Richtmyer  (Ref.  35,  p,  220)  for  this  case. 
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In  the  case  |i  =  0  it  is  also  possible  to  derive  a  stability  criterion  from  (B-19) 
and  (B-21)  under  the  condition  that  l-  a  Ax,  where  a  is  some  constant,  instead 
of  holding  l  constant.  Then  the  restriction  is  found  to  be  relaxed  to 


(B-24) 


Criterion  (B-23)  holds  in  the  region  of  the  shock,  and  criterion  (B-24)  is  valid 
in  regions  away  from  the  shock.  However,  if  true  viscosity  is  present  (u  ^  0) 
then  we  must  always  have  At/(Ax)  ■  0(1)  for  stability;  if  At/ Ax  ■  0(1)  one  of 
the  eigenvalues  goes  to  infinity.  This  results  since  in  this  case  the  term 
in  (B-19)  does  not  have  a  factor  of  order  (Ax)^  but  of  order  unity.  ° 
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Figure  1  Crater  profile  vs  impact  velocity  for  lead  target  and  various 
projectile  materials. 
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igure  t  Crater  profile  vs  impact  velocity  for  copper  target  and  various 
projectile  materials. 
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Figure  3  Crater  profile  vs  impact  velocity  for  steel  target  and  various 
projectile  materials. 
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Figure  4  Crater  profile  vs  impact  velocity  for  aluminum  target  and 
various  projectile  materials. 


Figure  5  Penetration  parameter  vs  impact  velocity  for  steel  target  and 
projectile. 
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Figure  7 


Penetration  parameter  vs  impact  velocity  for  lead  target  and 
projectile. 
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Figure  9 


Crater  volume  parameter  vs  impact  velocity  for  steel  target 
and  projectile. 


(Vc/Vp) 


VELOCITY  (  cm/a  sec.) 


Figure  10  Crater  volume  parameter  vs  impact  velocity  for  aluminum 
alloy  target  and  projectile. 
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Figure  11  Crater  volume  parameter  v8  impact  velocity  for  lead  target 
projectile. 
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Figure  12  Crater  volume  parameter  vs  impact  velocity  for  copper  target 
projectile. 
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(b) 


Figure  13.  Schematic  representations  of  the  (a)  forces  in  a  Bingham  body  and  (b) 
the  dependence  of  the  viscosity  on  the  strain-rate. 
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«  STRAIN 


Figure  14  Illustration  of  the  quantities  in  Malvern's  strain-rate  dependent 
constitutive  relation. 
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Figure  15  Illustration  of  impact  situation  in  (a)  laboratory  coordinates  and 
(b)  center  of  mass  coordinates. 
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Figure  17  Display  of  (a)  impact  situation  immediately  after  impact  and  (b) 
the  nomenclature  used  to  describe  the  initial  and  boundary  data. 
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Figure  18  Rankine-Hugoniot  pressures  calculated  for  two  semi-infinite 
bodies  of  indicated  material  impacting  at  velocity  v  . 
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Figure  19  Rankine-Hugoniot  internal  energies  calculated  for  two  semi- 
infinite  bodies  of  indicated  material  impacting  at  velocity  vQ. 
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Figure  20  Rankine-Hugoniot  densities  calculated  for  two  semi -infinite 
bodies  of  indicated  material  impacting  at  velocity  v  . 
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Figure  21  Hugoniot  curves  for  indicated  material. 
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Figure  22  Display  of  (a)  Eulerian  space-mesh  points  and  (b)  initial  and 
boundary  data  for  the  implicit  difference  scheme. 
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Figure  24.  Calculated  pressure  profiles  compared  with  the  Rankine-Hugoniot  sol 
tion(Casel:  v  =  .5,  t  =0). 
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Figure  25.  Pressure  profiles  (Case  2:  v  ■  .5,  t  *  0,  u  ■  08) 
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Figure 26.  Pre 


Figure  26.  Pressure  profiles  (Case 


Figure  27. 
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Figure  28.  Pressure  profiles  (Case  7:  v  ■  .5,  To  -  .01,  ■  .8). 
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Figure  29.  Pressure  profiles  (Case  11:  v  ■  .5,  t  -  .  1,  (i  =  .8). 
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igure  32  Pressure  profiles  (Case  32: 
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Figure  33  Pressure  profiles  (Case  48:  vQ  =7.5,  =1,  =  8), 
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PRESSURE (MEGABARS) 
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Figure  35.  Pressure  profiles  (Case  52:  vq  ■  1,  tq  =  0,  mq  *  8). 
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VON  MISES  FLOW  STATISTIC  { MEGABAR 
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Figure  56.  Von  M isos  flow  statistic  profiles  and  the  envelopes  of  their  maximum  values. 
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Figure  37  Thickness  of  the  stable  profile  as  a  function  of  impact  velocity. 
The  effect  of  the  strength  term  is  much  less  than  that  of  the 
viscosity. 


99 


IM  PACT  VELOCITY  (  CM  /  fi  SE 


Velocity  of  the  stable  pressure  profile  as 
impact  velocity.  The  viscosity  has  practi 


Figure  39  The  pressure  behind  the  stable  pressure  profile  as  a  function  of 
impact  velocity.  The  viscosity  has  very  small  effect. 
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Figure  42.  Profiles  of  stress  normal  to  the  wave  front  (Case  15:  vq  =  .  5,  tq  =  1,  Hq=.  8). 
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